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Abstract 

Response  of  a  shock  front  to  small  pre-shock  non-uniformities  of  density, 
pressure  and  velocity  is  studied  theoretically  and  numerically.  These  pre-shock  non¬ 
uniformities  emulate  imperfections  of  a  laser  target,  due  either  to  its  manufacturing,  like 
joints  or  feeding  tubes,  or  to  pre-shock  perturbation  seeding/growth,  as  well  as  density 
fluctuations  in  foam  targets,  “thermal  layers”  near  heated  surfaces,  etc.  Similarly  to  the 
shock-wave  interaction  with  a  small  non-uniformity  localized  at  a  material  interface, 
which  triggers  a  classical  Richtmyer-Meshkov  (RM)  instability,  interaction  of  a  shock 
wave  with  periodic  or  localized  pre-shock  perturbations  distributed  in  the  volume  distorts 
the  shape  of  the  shock  front  and  can  cause  a  RM-type  instability  growth.  Explicit 
asymptotic  formulae  describing  distortion  of  the  shock  front  and  the  rate  of  RM-type 
growth  are  presented.  These  formulae  are  favorably  compared  both  to  the  exact  solutions 
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I.  INTRODUCTION 


When  an  initially  symmetric  shock  wave,  perfectly  planar  or  spherical,  propagates 
into  a  non-uniform  medium,  its  front  shape  becomes  distorted.  The  perturbed  shock  front 
leaves  a  sonic/entropy/vorticity  perturbation  field  in  the  shocked  fluid  behind  it. 
Perturbations  produced  this  way,  in  turn,  can  trigger  either  a  Rayleigh-Taylor  (RT)  or  a 
Richtmyer-Meshkov-like  (RM-type)1  instability  growth,  particularly  if  an  ablation  front 
or  a  material  interface  is  present  nearby. 

The  response  of  shock  fronts  and  post-shock  flows  to  the  pre-shock  non¬ 
uniformities  is  of  interest  for  a  wide  range  of  applications,  primarily  in  the  area  of  laser 
fusion: 

1)  In  direct-drive  laser  fusion,  one  of  the  most  promising  methods  of  suppressing 
the  laser  imprint  and  mitigating  the  RT  growth  through  tailoring  the  target  adiabat  is 
based  on  irradiating  the  target  with  a  short  laser  pulse  (spike)  prior  to  the  main  pulse  that 
accelerates  and  implodes  it.  The  spike  produces  a  decelerating  shock  wave  propagating 
into  the  target,  and  a  rarefaction  wave  immediately  following  it.  Such  “impulsive 
loading”  of  the  target  has  been  shown  to  produce  strong  oscillations  of  the  shock  front 
and  the  areal  mass  in  the  shocked  material,  amplifying  the  non-uniformities  due  to  the 
target  manufacturing  imperfections  and  to  the  spike  pre-pulse  imprint.4’ 3  The  subsequent 
strong  shock  wave  produced  by  the  main  laser  pulse  propagates  through  the  perturbation 
field  left  after  the  spike. 

2)  In  some  of  the  direct-drive  laser  fusion  target  designs,  the  ablator  is  a  layer  of 
DT  wetted  or  empty  plastic  foam6  that  helps  increase  the  driving  laser  beam  absorption  in 
the  ablator  and  improve  its  stability.  Foam,  being  uniform  on  the  average,  contains  small- 
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scale  density  non-unifonnities,  which  might  affect  the  macroscopic  shock  dynamics.7  It  is 
important  to  know  whether  or  not  the  small-scale  foam  non-unifonnities  can  also 
contribute  to  the  shock  and  areal  mass  perturbation  growth  in  the  target. 

3)  There  are  certain  practical  advantages  of  having  spherical  laser  fusion  targets 
manufactured  with  some  localized  non-uniformities,  such  as  a  “thin  annular  joint 
extending  from  the  inner  to  the  outer  shell  radius,”  where  two  beryllium  hemispheres  are 

o 

bonded,  or  “a  narrow  bore  fill  tube  to  charge  an  ignition  capsule  in  situ  with  deuterium- 
tritium  fuel.”9 

4)  Deformation  of  a  blast  wave  in  long  thin  regions  of  lower  density  was  first 
observed  during  nuclear  testing  and  later  labeled  “the  thermal  layer  effect,”10  and 
extensively  studied,  see  Ref.  1 1  and  references  therein.  Such  pre-shock  “thermal  layers,” 
where  the  temperature  is  higher  and  the  density  lower  than  in  the  ambient  medium, 
thereby  forming  a  “wall-supported  shock  precursor,”  can  be  produced  near  heated  walls 
in  shock-tube,12  laser13  and  plasma-focus14  laboratory  experiments. 

Small-amplitude  theory  describing  the  interaction  of  a  shock  wave  with  a  non¬ 
uniformity  localized  at  a  material  interface  had  been  first  developed  by  Richtmyer15  and 
later  advanced  by  many  researchers,  see  Refs.  16-17  and  references  therein.  Given  the 
importance  of  shock  interaction  with  small  pre-shock  non-uniformities  for  many 
applications,  it  is  rather  surprising  that  Richtmyer’  approach  has  not  been  extended  so  far 
to  produce  a  theory  describing  the  shock  distortion  and  the  post-shock  flow  evolution 
caused  by  the  pre-shock  non-uniformities  distributed  in  the  volume. 

In  this  article,  we  present  such  a  theory  based  on  small-amplitude,  linear  stability 
analysis,  for  a  planar  geometry.  The  background  post-shock  and  pre-shock  fluids  are 
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assumed  to  be  uniform,  and  the  unperturbed  shock  strength  maintained  constant  in  time. 
The  pre-shock  perturbations  considered  are  non-propagating  structures,  such  as 
density/entropy  modulation,  as  well  as  standing  sonic  waves.  We  present  the  results  for  a 
shock  wave  which  is  either  incident  from  a  half-space  where  there  is  no  non-uniformity 
(isolated)  or  supported  by  a  piston.  A  free  surface  and  an  ablation  front,  where  a  constant 
pressure  is  maintained,  and  a  rigid  wall  are  examples  of  the  piston  considered  below. 
Some  of  the  results  of  our  theory  are  directly  applicable  to,  and  the  other  are  easily 
generalized  for  the  case  when  a  RM-unstable  material  interface  serves  as  a  piston  driving 
a  constant-strength  transmitted  shock  wave  into  a  weakly  non-unifonn  fluid.  Our  theory 
applies  to  the  shock  interaction  with  small-amplitude  non-uniformities,  both  periodic  and 
localized  in  space,  such  as  “joints”,  “tubes”  and  “thermal  layers.”  It  is  not  applicable  to  a 
blast  wave,  whose  strength  decreases  with  time,11  or  to  a  shock  wave  driven  up  the 
density  gradient  produced  by  a  spike.”'  Further  development  of  the  theory  to  cover  these 
cases  is  obviously  needed,  which  is  beyond  the  scope  of  the  present  article. 

In  Section  II  we  present  formulation  of  the  problem,  perturbation  equations, 
boundary  conditions,  and  the  elementary  theory  of  shock  distortion  for  the  simplest  cases. 
In  Section  III  we  describe  the  general  theory  and  compare  the  time  evolution  of  the 
perturbed  shock  front  to  our  explicit  asymptotic  fonnulae.  Section  IV  compares  our 
analytical  results  with  numerical  simulations.  In  Section  V  we  conclude  with  a 
discussion. 
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II.  FORMULATION  OF  THE  PROBLEM  AND  ELEMENTARY  THEORY 


A.  Perturbation  equations,  initial  and  boundary  conditions 

Consider  a  planar  shock  wave  incident  at  t  =  0  upon  a  half-space  x  >  0 .  The  fluid 

ahead  of  the  shock  front  is  weakly  non-uniform: 

P  =  Pi  +5p(x,y,0;  p  =  px+bp(x,y,t); 
vv  =  8vv  (x,  y,  t);  vy  =  Svv  (x,  y,  t); 

we  suppose  that  Sp  /  p , ,  Sp  /  p] ,  8vx  /  ax  and  8v  / ax  are  all  much  less  than  unity,  where 


ax  is  the  pre-shock  speed  of  sound  (indices  1  and  2  refer  to  pre-  and  post-shock  quantities, 
respectively).  In  the  zero-order  approximation,  neglecting  the  small  non-unifonnities,  we 
have  a  planar  shock  wave  propagating  at  constant  velocity  D  into  a  uniform  fluid.  The 
mass  velocity  of  the  fluid  behind  the  shock  front  with  respect  to  the  pre-shock  fluid  is 
denoted  U;  obviously,  the  density  compression  ratio  in  the  shock  wave  is  related  to  these 
velocities  by 


R-p2_  D 
Pi  D-U' 


(2) 


Below  we  assume  ideal  gas  equation  of  state  of  the  fluid  with  adiabatic  exponent  y.  The 

1  R 

average  post-shock  values  of  density,  pressure  and  sound  speed  are 


P  2 


a  2 


(y  +  l)Mf 

{y-\)M2x+  2Pl’ 


2y M:  -y  +1 

Pi  = - - - Pi, 

y  +  1 


^(2yM12-y  +  l)[(y-l)M12+2] 
(y+l)Mi 


(3) 


where  Mx  =  D!  ax  >  1  is  the  Mach  number  of  the  incident  shock.  We  will  use  the 
reference  frame  in  which  the  shocked  fluid  is  at  rest  and  x  =  0  corresponds  to  the 
position  of  the  planar  shock  front  at  t  =  0  .  In  this  reference  frame,  the  shock  front 
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propagates  at  the  velocity  D  —  U  with  respect  to  the  shocked  fluid.  The  corresponding 
downstream  Mach  number  is 


M, 


D-U  _  (y-l)M12+2  <  ] 


2yMf  -y  +  1 


(4) 


Our  small-amplitude  analysis  is  based  on  Fourier  decomposition  of  all 
perturbations.  We  assume  their  dependence  on  the  transverse  coordinate  y  in  the  form 
exp(z'ky)  ,  where  k  =  k  =  2%  / Xv ,  see  Fig.  l(a)-(c).  Below  k  without  subscript  always 
means  k  ,  and  the  common  factor  exp  {iky)  in  perturbation  variables  is  omitted.  Let  us 


introduce  dimensionless  first-order  perturbation  variables 
~  _  Sp  _ _  8/?  _  1  6/7  _  __  Sv\  _  .  6vv 


p2^2  Y  Pi 


vy=l~ 


(5) 


Then  the  linearized  ideal  fluid  equations  are  written  as 
1  dp  dv 

ys7+a7+fa,'=0. 


(6) 


1  dp  dv  . 

-^  +  -±  +  kv=0, 


a2  dt  dx 


(7) 


1  5vv.  +  ^P=  o 


a2  dt  dx 


(8) 


1  _ 

-kp  =  0 . 


a2  dt 


(9) 


Excluding  velocity  perturbations  from  (6)-(9),  we  obtain  the  Klein-Gordon 
equation  for  the  dimensionless  pressure  perturbation: 

1  d2p  d2p  :  ,  2 


a2  dt 2  dx2 


+  k~  p  =  0 . 


(10) 
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Equations  (6)-(10)  must  be  supplemented  by  the  boundary  conditions  at  the 
surface  of  the  piston  driving  the  shock  and  at  the  shock  front.  For  both  the  rigid  piston 
and  the  free  surface  boundary  conditions,  the  piston  is  formed  by  the  fluid  particles  which 


initially  were  located  at  x  =  0 .  For  the  rigid  piston  the  boundary  condition  at  x  =  0  is 
vr=0;  (11) 


for  the  free  surface,  it  is 

p  -  0 . 


(12) 


If  the  shock  wave  is  driven  by  a  laser,  as  in  laser-fusion-relevant  conditions,  then 
the  boundary  condition  should  be  formulated  at  the  time-dependent  position  of  the 
ablation  front,  at  x  =  vat ,  where  the  ablation  velocity  va  =  ml  p2 ,  and  m  is  the  mass 

ablation  rate.  The  corresponding  Mach  number  is  defined  as  M a  =  va  la2 ;  in  most  cases 
it  is  small,  M a  «  1 ,  because  the  plasma  flow  near  the  ablation  front  is  subsonic,  almost 

isobaric.  Below  we  approximate  the  ablation  front  with  a  sharp  boundary  and  use  the 
boundary  conditions  due  to  Sanz  and  Piriz,19  formulated  as  in  Eqs.  (B7)  and  (9)  of  Ref. 
20: 


P  +  ^  ^  M“  khxa  +  Ma  [(1  -  rD)vx  +  +  =  0  .  (13) 

rD  ka2  at 

— ^5x„,  +  Makhxa  -  vx  =  0 .  (14) 

a2  dt 

Here,  8xa  is  the  ripple  amplitude  of  the  ablation  front  and  rn  =  pw  /  p2  is  the  blowoff- 
plasma-to-ablation  front  density  ratio,  which  also  is  typically  small,  rD«  1 .  In  the  limit 
Ma  — >  0 ,  rD  — >  0  ,  M 1 1  rD  — >  0  Eqs.  (13)  and  (14)  are  reduced,  respectively,  to  the  free- 
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surface  boundary  condition  (12)  and  to  the  conventional  relation  between  the  normal 
velocity  and  displacement  amplitudes  at  the  free  surface. 

Finally,  the  “isolated  shock”  assumption  means  that  the  incident  shock  wave 
arrives  to  the  x  =  0  boundary  from  a  half-space  filled  with  a  uniform  fluid.  A 
weak/contact  discontinuity  is  instantly  formed  at  the  fluid  particles,  where  the  density 
modulation  begins  (which  correspond  to  x  =  0  in  the  shocked  frame  of  reference).  The 
perturbations  of  pressure  p  and  the  longitudinal  velocity  vx  are  continuous  at  x  =  0 , 
whereas  the  relative  density  perturbation  p  and  the  x-derivative  of  the  lateral  velocity, 
vv ,  undergo  a  discontinuity.  The  flow  in  the  reflected  sonic  wave  is  governed  by  the 

same  equations  (6)-(9),  which  need  to  be  solved  additionally  to  the  same  equations 
describing  the  fluid  between  x  =  0  and  the  shock  front  because  of  the  discontinuity 
at  x  =  0 .  The  boundary  conditions  at  the  planar  front  of  the  reflected  sonic  wave,  at 
x  =  —a2t ,  are 

vx+p  =  o,  v  v  =  0  .  (15) 

The  boundary  conditions  at  the  shock  front  are  perturbed  Rankine-Hugoniot 
(RH)  conditions  derived  in  Ref.  18,  §90.  We  write  them  down  here  for  a  particular  case 
of  density/entropy  pre-shock  modulation  at  constant  pressure  and  zero  velocity, 
characterized  by  the  amplitude  8pA(x)  ,  which  is  supposed  to  be  relatively  small  at  all  x: 


|8pi.(x)|/p1  «1. 


A  -M2(R-  l)k8xs  =  0  . 


27 


(Mi+1L_  n5P* 


M~M2 


-p  =  -M2(R-l)- 


(16) 

(17) 
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(18) 


P  - 


sp* 


M;Ml '  Pl 


2  d  5 

- ox 

a2  dt 


(Y  +  1)^_  ^  d5P* 


2  M, 


-p  =  —M2R- 


(19) 


Here,  8xs  is  the  displacement  amplitude  of  the  shock  front,  R  and  M2  are  given 
by  (2)-(4).  In  the  reference  frame  where  the  shocked  fluid  is  at  rest,  the  boundary 
conditions  (16)-(19)  are  imposed  at  xs  =  (D-U)t  =  M2a2t ,  which  corresponds  to 
x  =  Dt  in  the  reference  frame  where  the  pre-shock  fluid  is  at  rest.  The  initial  coordinate- 
dependent  pre-shock  density  modulation  profde  8pt(jc)  thereby  translates  into  a  time- 

dependent  term  8pt(*)/  Pi  in  (17)-(19).  Without  this  term,  the  boundary  conditions  (17)- 

(19)  are  the  same  as  used  by  all  the  researchers  who  studied  propagation  of  rippled  shock 
waves  into  a  uniform  ideal  gas,  starting  from  Refs.  21-22.  If  the  pre-shock  fluid  also 
contains  perturbation  of  pressure  and/or  velocity,  as,  for  example,  sonic  waves  (see 
Section  III.D),  the  corresponding  driving  terms  should  be  added  to  the  right-hand  sides  of 
Eqs.  ( 16)-(  19). 

Figure  1  demonstrates  some  examples  of  pre-shock  density  modulation  profiles 
studied  in  this  article.  The  simplest  case  is  periodic  lateral  variation,  independent  from 
the  longitudinal  coordinate  x,  “stripes,”  Fig.  1(a),  8p  °c  cos(kvy)  ,  which  corresponds  to 

constant  values  of  the  density  modulation  terms  in  the  right-hand  sides  of  the  boundary 
conditions  (17)-(19).  Periodicity  in  both  lateral  and  longitudinal  directions,  a  “quilt”  of 
Fig.  1(b),  8p  oc  cos(k  v)cos(kvx),  corresponds  to  single-mode  oscillations  of  these  terms 


in  ( 1 7)-(  1 9)  with  the  frequency  co  =  kxD  .  Figure  1(c)  demonstrates  the  case  where 
density  modulation  is  periodic  in  the  lateral  direction  but  concentrated  near  the  x  =  0 
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surface,  8p  °<  cos(Ay)exp(-gx)  ;  the  corresponding  density  modulation  terms  in  Eqs. 

( 1 7)-(  1 9)  decay  with  time  as  cxp(-qDl) .  Finally,  Figure  1(d)  shows  a  Gaussian  density 
profile  Sp  °c  exp(-g2 y2) ,  emulating  a  “thermal  layer”  or  a  “joint”.  This  localized  pre¬ 
shock  density  modulation  is  not  periodic  in  the  lateral  direction;  hence  our  single-mode 
equations  (6)-(  1 0),  (13)-(  19)  are  not  directly  applicable  to  study  it.  Nevertheless,  having 
solved  the  perturbation  problem  for  the  simplest  “stripes”  case  of  Fig.  1(a),  we  can 
address  the  ’’joint”  case  of  Fig.  1(d)  by  using  Fourier  integral,  in  the  same  way  as  done  in 
Refs.  23  for  a  shock  wave  driven  by,  or  reflected  from,  a  rigid  piston  with  a  Gaussian 
perturbation. 

Figure  2  illustrates  the  perturbation  problem  formulation  for  the  “stripes”  case  of 
Fig.  1(a)  and  the  shock  wave  driven  by  a  rigid  piston.  The  shock  front  moves  from  the 
piston  in  the  positive  x  direction.  In  our  reference  frame,  where  the  piston  is  at  rest,  the 
unperturbed  coordinate  of  the  shock  front  is  xs  =  (D-U)t  =  M2a2t ,  see  above.  Deviation 
of  the  shock  front  from  its  original  planar  shape  is  caused  by  the  interplay  of  the  pre¬ 
shock  non-uniformity  of  density  and  the  post-shock  non-uniformity  of  all  the  fluid 
variables.  Our  goal  is  to  determine  the  time  dependence  of  the  shock  front  displacement 
amplitude,  Sxv  (/) ,  as  well  as  the  dynamics  of  the  other  perturbed  variables  in  the  shocked 
fluid. 

B.  Elementary  theory 

The  problem  of  shock  interaction  with  pre-shock  density  perturbation  has  an 
elementary  asymptotic  solution  for  the  simple  case  of  an  isolated  shock  wave  coming 
from  the  unperturbed  fluid.  In  this  case,  it  is  sufficient  to  calculate  the  amplitudes  of  the 
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perturbations  which  represent  waves  radiated  downstream  with  respect  to  the  shock  front. 
The  solution  is  obtained  for  the  general  case  of  a  “quilt,”  Fig.  1(b): 


— -  =  £ k  cos (kxx)  ,  (20) 

Pi 

where  ek  «  1  is  a  dimensionless  constant.  The  “stripes”  case  of  Fig.  1(a)  is  recovered  in 
the  limit  kx  — »  0  . 

Assuming  that  the  displacement  of  and  the  pressure  perturbation  at  the  shock  front 
after  some  transient  period  oscillate  harmonically  at  the  frequency  co  =  kxD  ,  we  can  find 

the  amplitudes  of  these  oscillations.  To  do  this,  we  decompose  the  perturbation  behind 
the  shock  front  into  an  outgoing  sonic  wave,  a  vortex  wave  and  an  entropy  wave: 


p=r\  v,=v»+vw,  y,=v«+v«,  p=p«+pM 


where  the  superscripts  5,  v  and  e  refer  to  the  sonic,  vorticity  and  entropy  perturbations, 
respectively. 


The  sonic  pressure  perturbation  p  =  /A'1  satisfies  the  wave  equation  (10).  We  seek 


a  solution  of  (10)  that  at  x  =  M2a2t  oscillates  at  the  frequency  to  in  the  form 


p{s)  =  p0  exp[-/co?  +  kr|(x-M2<72/)], 


where  p  is  a  constant  to  be  determined.  Substituting  (22)  into  (10),  we  obtain  the 
dispersion  relation 

(/co  +  ka2M2r\  )2  - k2ajr\ 2  +  k2a \  =  0 , 


from  which  we  find  the  parameter  T): 


J= (±)n 


1  -Mi 


0  0  )  ’ 
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where  we  have  defined 


c 


0  — 


kxRM  2 

ky^-M22 


(25) 


We  see  that  the  sonic  perturbations  downstream  behave  differently  for  the  values 
of  the  parameter  C,  0  less  and  greater  than  unity.  In  the  former  case,  the  real  part  of  the 

parameter  T|  is  nonzero,  which  implies  that  the  corresponding  pressure  perturbations  (22) 
either  grow  or  decay  exponentially  with  the  distance  from  the  shock  front.  Physically,  it 
means  that  our  combination  of  the  frequency  co  =  k  v  D  and  the  lateral  wavenumber  kv 

does  not  correspond  to  a  running  sonic  wave.  The  perturbation  amplitudes  are  thereby 
maintained  finite  at  the  shock  front  by  the  external  driving  force  (pre-shock  density 
modulation)  and  decay  with  the  distance  from  the  shock  front.  The  corresponding 
condition  on  the  ratio  of  longitudinal  and  lateral  wavenumbers  can  be  re-written  as 


k 

kv  RM2 


(26) 


If  the  long-wavelength  condition  (26)  is  satisfied,  and  thereby  C,  „  <  1  >  then  we  are 
interested  in  the  sonic  mode  that  decays  with  the  distance  from  the  shock  front,  at 
x  — >  — °o  ,  that  is,  one  corresponding  to  the  solution  (24)  with  the  positive  real  part 

J—  (vnfi+iM2e„).  <27) 

^-M2 


If  (26)  is  not  satisfied,  C,  0  >  1  >  then  both  solutions  (24)  are  imaginary,  and  the 
pressure  perturbations  are  running  sonic  waves.  The  outgoing  wave  corresponds  to 


(-V^M+m2c0). 


(28) 
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Sonic  velocity  perturbations  vx(,s)  and  v^s)  are  proportional  to  the  same  exponential 


factor  (22),  and  their  amplitudes  are  found  from  Eqs.  (8),  (9): 


y(s) 

Vo 


Vs) 


P  o 


M2r\+it;J\-M22  ' 


(29) 


Velocity  perturbations  in  the  vortex  wave  are  arbitrary  functions  of  x,  satisfying 
the  divergence-free  condition  3v)(v)  / c)x  +  kv(}v>  =  0 .  At  the  shock  front,  x  =  M2a2t ,  these 

functions  should  oscillate,  which  detennines  the  coordinate  dependence  of  these 
velocities.  For  the  longitudinal  velocity  we  have 


^v(v)  =v,(o)exp 


cox 


-7 - 

Y  M2a2  j 


~( v) 

VoexP 


-ikx 


M, 


(30) 


The  coordinate  dependence  is  the  same  for  the  lateral  velocity  perturbation,  whereas  its 


amplitude  vj,  is  expressed  via  vxy  with  the  aid  of  the  divergence-free  condition: 


V (v)  =  7 
Vy  0  1 


M, 


2_v(v) 

yx0  • 


(31) 


At  the  shock  front  p  =  p0 ,  vx  =  v),)'  +  ,  v  =  .  The  two  unknown 


amplitudes,  p0  and  vy  ,  and  the  dimensionless  shock  ripple  amplitude  kSxs  are 
immediately  detennined  from  the  RH  boundary  conditions  (16),  (17)  and  (19),  where  we 
have  to  substitute  1  /(ka2)(d  / dt  )  =  -i%o2Jl-M2  . 

For  the  long-wavelength  “quilt”  case,  when  (26)  is  satisfied,  the  asymptotic 
expression  for  the  shock  ripple  amplitude  is 


k8xs(t)  =  -£kQj  1  Pl?0" 


l+P2Co2  +  P3Cc 


cos(  k  x.Dt  —(f)). 


(32) 
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where  the  constants  are  defined  as 


Q 


P2 


RM : 


P.=l- 


M\M\  ’ 


f 


2M~  —  1  — - 


V 


M, 


P3 


1  +  - 


M: 


-4  M 


2  9 


(33) 


2  m2mX0Ji^ 

<|)  =  arctan — — — p  "  -  arctan 


2 
5  0 


M?-(M;+  1)C02 


and  the  value  of  arctan(...)  in  the  definition  of  (f)  is  taken  between  n  12  and  7t  if  its 
argument  is  negative. 

In  particular,  for  the  limiting  case  of  “stripes”  pre-shock  density  modulation, 
kx  =  0 ,  hence  £  0  =  0  ,  and  (32)  reduces  to 


kbxs(t)  =  -zkQ.  (34) 

Physical  interpretation  of  (34)  is  simple.  When  a  planar  shock  wave  starts 
propagating  along  the  “stripes”  of  the  density  profile  shown  in  Fig.  1(a),  initially  it  moves 
faster  where  the  density  is  lower.  It  happens  because  the  post-shock  flow  is  subsonic,  and 
the  pressure  there  rapidly  equilibrates,  and  the  shock  velocity  scales  as  D  °c  ^fpTp  , 
where  p  is  the  post-shock  pressure,  and  p  is  the  pre-shock  density.  The  initially  planar 
shape  of  the  shock  front  is  distorted,  becoming  convex  where  the  density  is  lower  and 
concave  where  it  is  higher,  see  Fig.  2,  and  thereby  equilibrating  the  pressure  and  reducing 
the  difference  in  velocities.  Eventually  the  shock  approaches  a  constant  velocity  and  a 
steady-state  rippled  shape,  the  asymptotic  shock  ripple  amplitude  being  given  by  (37). 
Since  Q  >  0 ,  the  negative  sign  in  (34)  means  that  the  shock  is  farther  from  the  piston 
where  the  pre-shock  density  is  initially  lower. 
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For  the  “quilts”  with  very  long  longitudinal  wavelengths,  kx  «  k  ,  C,0  «  1 ,  Eq. 
(32)  demonstrates  that  the  shock  ripple  amplitude  adjusts  to  the  slowly  varying  local 
value  of  the  density  modulation  amplitude,  £k  cos{kxDt).  The  coefficient  of 
proportionality  between  the  shock  ripple  amplitude  and  the  pre-shock  density  modulation 
amplitude  is  close  to  the  value  -  Q  characteristic  of  the  “stripes”  case  (the  correction  is  of 

the  order  of  ^02).  There  is,  however,  a  phase  delay  (])  between  the  pre-shock  density 
modulation  and  the  shock  ripple,  so  that  the  former  peaks  not  when  the  shock  passes  a 
valley  of  density  modulation,  but  slightly  later.  With  increased  kx ,  the  amplitude  of 
oscillations  monotonically  increases  compared  to  the  “stripes”  case,  by  a  factor  varying 
from  1  at  ^  0  — >  0  to  1  /  M2  >  1  at  £  0  =  1 .  The  phase  shift  (])  also  increases  from  0  to  n  12 


at  C0  =  1- 


For  the  short- wavelength  “quilt,”  when  (26)  is  not  satisfied  and  C,  0  >  1 ,  the 


asymptotic  expression  for  the  shock  ripple  amplitude  is 

C„ 


kbxs(t)  =  ~zkQ 


M, 


-1  -M: 


-sin  (kxDt).  (35) 


In  this  regime,  the  phase  shift  between  the  density  modulation  and  the  shock  ripple  has  a 
constant  value  of  7t  /2  ,  so  that  the  shock  ripple  amplitude  peaks  where  the  pre-shock 
density  modulation  has  zero  value.  The  asymptotic  normalized  amplitude  (35)  of  the 
short-wave  shock  ripple  oscillations  decreases  with  increased  longitudinal  wavenumber 
as  k  /  kx,  and  the  limit  kx  /  ky  — >  °°  it  can  be  presented  as 


khxs(t)  =  -£k 


M\M2  + 1 


2Mj  M2  +m;  +1 


sin  (kxDt). 


(36) 
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The  above  formulae  have  been  obtained  for  the  isolated  shock  boundary 
condition.  However,  for  the  “stripes”  and  long-wavelength  “quilt”  cases  the  outgoing 
sonic  perturbations  exponentially  decay  with  the  distance  from  the  shock  front.  Therefore 
the  perturbed  shock  front  in  these  cases  becomes  decoupled  from  the  piston  at  late  time, 
t  —>  oo .  Consequently,  the  asymptotic  formulae  (32)  and  (34)  are  valid  for  arbitrary 
boundary  conditions  at  the  piston. 

This  is  not  the  case  for  the  short-wavelength  regime,  when  the  perturbed  shock 
front  radiates  running  sonic  waves  into  the  shocked  fluid.  If  the  incident  shock  wave  is 
not  isolated  but  driven  by  a  piston,  then  these  sonic  waves,  having  been  reflected  from 
the  piston,  might  reach  the  shock  front  again,  adding  new,  Doppler-downshifted 
frequencies  to  its  oscillations,  and  thereby  making  the  asymptotic  formulae  (35)  and  (36) 
inapplicable.  Extending  the  above  derivation  to  cover  this  case  is  possible  but  probably 
too  complicated.  To  study  the  short-wavelength  case  in  the  presence  of  a  piston,  we  will 
apply  more  advanced  theoretical  methods  described  in  Section  III. 

Here  we  limit  ourselves  to  the  simplest  case  of  longitudinal  density  stratification 
that  produces  sonic  waves  propagating  from  the  piston  and  back.  Let  us  take  the  lateral 
wavenumber  k  =  0  and  Sp  /  =e0/(jc),  where  £0  «  1  is  the  small  dimensionless 

amplitude  of  the  pre-shock  longitudinal  density  modulation,  and  / (x)  is  an  arbitrary 
function  describing  its  profile.  Then  Eq.  (10)  becomes  a  one-dimensional  wave  equation, 
whose  general  solution  is  a  combination  of  running  waves.  We  present  it  in  the  form 

v  =  —F(aJ  +  x)  +  RnF(aJ  -  x),  v  =  0, 

^  p  ^  y 

p  =  F(a2t  +  x)  +  RpF(a2t  -  x) ,  p  =  p  +  G(x) . 
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Here,  R  is  the  reflection  coefficient  for  the  sonic  waves  incident  upon  the  piston; 


R  =  1 ,  - 1 ,  and  0  for  the  rigid  piston,  free  surface  and  isolated  shock  boundary 

conditions,  (1 1),  (12),  and  (15),  respectively.  Obviously,  the  solution  (37)  identically 
satisfies  Eqs.  (6)-(  1 0)  and  the  boundary  conditions  (1 1),  (12)  or  (15)  with  arbitrary 
functions  F(x )  and  G(x) .  These  functions  are  easily  detennined  if  we  substitute  (37) 
into  (17)  and  (18).  The  result  is: 

F(x)  =  (R  -  1)(R,  + 1)£  (*,*,  y  f 

4  11=0 

where 


RM-,x 


1  +  M, 


1  -m2 

1  +  M, 


\«  A 


(38) 


R  = 


1  -m2 

1  +  M, 


2(y-l)M2-3  +  y 


(39) 


2(y  -1)M2  +3-y 

v  / 

is  the  reflection  coefficient  for  the  sonic  waves  normally  incident  upon  the  shock  front 
from  downstream.  Substituting  the  solution  (38)  into  the  boundary  condition  (18),  one 
can  easily  determine  the  entropy  perturbation  profile  G(x )  .  The  infinite  series  (38) 
describes  multiple  reverberations  of  sonic  waves  between  the  shock  front  and  the  piston; 
the  coefficient  [(1  -  M2 )  /(I  +  M2 )]"  accounts  for  the  Doppler  shift  of  the  reverberating 
wave  after  n  reflections  from  both  surfaces.  The  reflection  coefficient  R,  varies  between 
zero,  for  weak  shock  waves,  and  some  peak  value,  which  for  y  <  5/3  is  achieved  for 
strong  shock  waves.  This  peak  value  is  small  for  typical  gammas,  varying  from  -0.0557 
for  y  =  5/3  to  -0. 139  for  y  =  7/5  .  This  is  why  the  series  (38)  rapidly  converge  for  not- 


too-small  gammas. 
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To  give  an  example,  consider  a  periodic  pre-shock  density  stratification, 
f(x )  =£0  cosk.x  .  Substituting  the  solution  (37)  and  (38)  into  (19),  we  find  the  time 
dependence  of  the  shock  displacement  amplitude: 


kx8xs(t)  =  — £p - ^ -  + 1 - sin krDt 


2M[M2  +M-  + 1 


+  eo^ 


{\  +  Rs)2{\  +  M2)2_  y/ 


8  M2R 


n—\ 


„  r,  1  +  ^7 


V 


\n- 1 

( 

f \-M2  ' 

n  ^ 

sin 

k Dt 

1 +  M, 

/ 

V 

f  2  7 

) 

(40) 


For  the  isolated  shock  cas qR  =  0  and  (40)  is  reduced  to  (36),  as  it  should  be. 


III.  SOLUTION  BY  MEANS  OF  LAPLACE  TRANSFORMS 

Here  we  present  general  equations  describing  the  temporal  evolution  of  the 
perturbations  in  the  space  filled  with  the  compressed  fluid  for  both  situations  of  long  and 
short  wavelength.  This  can  be  done,  in  principle,  for  arbitrary  boundary  conditions  at  the 
surface  x  =  0 .  The  Laplace  transform  technique  that  we  use  was  first  developed  and 
applied  to  the  study  of  perturbed  shock  evolution  by  Zaidel  in  Ref.  22.  Here,  as  had  been 
done  in  Refs.  16,  the  idea  is  to  derive  exact  analytical  expressions  for  the  Laplace 
transform  of  the  shock  front  pressure  perturbations  or  of  the  shock  front  ripple.  For  the 
isolated  shock  case,  which  represents  the  simplest  of  the  boundary  conditions  imposed 
downstream  of  the  shock  surface,  we  arrive  to  an  exact  closed  fonn  expression  for  the 
Laplace  transform  of  the  perturbations.  Inverting  this  Laplace  function,  we  can  follow  the 
evolution  of  almost  any  quantity  of  interest  in  the  compressed  fluid,  always  in  linear 
theory.  For  all  other  cases,  the  advantage  of  using  exact  expressions  for  the  Laplace 
transforms  of  the  perturbation  variables  lies  in  the  possibility  of  getting  accurate 
information  on  the  asymptotic  behavior  of  the  involved  functions.  The  exact  solutions 
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describing  the  shock  distortion  and  mass  flow  in  this  Section  I  are  obtained  using  either 
the  exact  inversion  of  the  Laplace  transform,  for  the  isolated  shock  case,  or  the  Taylor 
series  expansion  method  (see  Ref.  1  and  references  therein),  for  all  other  cases. 

A.  Pre-shock  density  modulation  periodic  in  the  lateral  direction  -  “stripes” 

We  start  with  the  simplest  of  the  boundary  conditions,  that  of  an  isolated  shock 
front.  That  is,  we  assume  that  there  is  no  piston  or  surface  reflecting  the  sound  waves 
back  to  the  shock  front  from  behind.  In  general,  the  shock  front  will  be  rippled  and  the 
corrugations  in  its  shape  will  induce  a  mass  flow  which  is  tangential  to  its  surface.  This 
lateral  mass  flow,  in  turn,  creates  a  perturbed  pressure  gradient  along  it.  As  a 
consequence,  the  shock  radiates  sonic  waves  downstream,  which  are  never  reflected  back 
to  the  shock  front.  This  case  represents  a  direct  extension  of  the  original  Richtmyer’s 
problem  (an  incident  planar  shock  interacts  with  an  interfacial  density  non-uniformity15) 
for  the  case  when  the  non-uniformity  is  distributed  in  the  volume.  In  laser-fusion-related 
applications,  the  shock  is  typically  driven  by  a  not-too-distant  piston.  Our  isolated  shock 
case  serves  as  a  basis  for  further  generalizations  that  take  into  account  the  boundary 
conditions  at  the  piston,  and  represents  a  good  approximation  in  the  limit  of  weak 
incident  shocks  for  arbitrary  boundary  conditions  at  the  piston. 

For  the  “stripes”  case  of  Fig.  1(a),  SpA,  /  pt  =  ek  =  const ,  where  ek  «  1  is  the 
small  amplitude  of  the  pre-shock  density  perturbation  exp  (iky)  .  As  in  Ref.  16,  we 

make  a  transfonnation  of  variables  from  x  and  t  to  r  and  0  ,  first  introduced  in  the 
context  of  a  similar  problem  by  Zaidel  in  Ref.  22  and  defined  by: 

r  coshQ  =  ka2t,  rsinh0=Ax.  (41) 
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The  dimensionless  shock  front  coordinate  0  v  is  expressed  via  the  downstream  shock 
Mach  number:  tanh0s  =  M2 .  It  is  convenient  to  work  with  the  Laplace  transform  with 

respect  to  the  variable  r.  We  define,  for  example,  the  Laplace  transform  of  the  shock 
pressure  perturbation  amplitude  as: 


P(s$s)  =  Ps(s)  =|  p(Qs,r)exp(-sr)dr,  (42) 

o 


where  the  variable  0  S  means  that  the  pressure  is  evaluated  at  the  shock  front  position. 

After  some  algebra,  which  can  be  adapted  to  this  problem  following  the  procedure 
outlined  in  Ref.  16,  see  Appendix,  we  arrive  to  an  exact  solution  for  the  Laplace 
transform  of  the  pressure  perturbations  at  the  shock  front  (assuming  that  there  are  no 
waves  reflected  back  to  the  shock  front,  as  also  discussed  in  Ref.  24)  in  the  “stripes”  case: 


p  ^  _  (l  oc10)s  ps()  +cc 


20 


(sV 5 2  +1  -a10s2  -an) 


(43) 


The  coefficients  a10,  an  and  a20  depend  on  the  fluid  equation  of  state  and  the  shock 
intensity.  The  coefficient  a20  here  is  different  from  that  defined  in  Ref.  16,  because  of  the 
different  nature  of  the  perturbations  at  the  shock  front.  In  Ref.  16,  the  shock  ripple 
amplitude  was  given  as  an  initial  condition  at  t  =  0+  ,  x  =  0 .  Here,  the  shock  front  ripple 
is  initially  zero,  and  becomes  nonzero  later  on,  due  to  the  density  modulation  in  the  pre¬ 
shock  fluid.  These  coefficients  for  an  ideal  gas  equation  of  state  (EOS)  valid  for  a  more 
general  case  of  a  “quilt”  pre-shock  density  modulation  are  given  by 


a 


(1  +  Mj2) 


10  2 M\M2 


a, 


2M-, 


a 


20 


■£kM  2 


y  + 1  2 


(44) 
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and  ps0  is  the  initial  pressure  perturbation  at  t  =  0 1  ,  at  the  surface  x  =  0  : 


e  M2(R-1) 

k  2(1 -a10)  ' 


(45) 


Here,  the  downstream  Mach  number  M2  and  the  shock  compression  ratio  R  are 
given  by  (4)  and  (2)-(3),  respectively,  and  C,0  -  by  (25);  in  our  “stripes”  case  C,0  =  0  . 

It  is  also  easy  to  obtain  the  solution  for  the  Laplace  transfonn  of  the  shock  front 
ripple  amplitude  k8xs  from  the  linearized  RH  equations  using  our  previous  results: 


kSXs(s)  = 
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4  Ml 
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+  -^-UJs  +1 -a 
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sV  S“  +  1— OC105  — OC  j  j 
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where  the  coefficient  a30  is  defined  as 
ekRM , 

«30  = - r - V- 

2  VI  ~Ml 


(47) 


From  Eqs.  (43)  and  (46)  we  can  obtain  the  temporal  evolution  of  the 
corresponding  perturbations.  A  conventional  way  of  doing  this  is  to  expand  the 
functions  given  by  Eqs.  (43)  and  (46)  as  a  series  in  the  variable  exp(-  sinlT1  s ).  This 
decomposition  for  the  Laplace  transform,  gives  rise  to  an  infinite  series  of  Bessel 
functions  in  the  time  domain.1617  24  However,  we  prefer  to  show  here  a  different 

approach,  that  of  explicitly  inverting  the  Laplace  functions  Ps  (5)  and  kSX  (s)  in  the 

complex  plane.  This  method,  albeit  apparently  more  complicated,  provides  us  with  an 
integral  representation  of  the  front  perturbations,  which  will  be  useful  for  many  other 
purposes,  as  for  example  the  estimation  of  the  late-time  asymptotic  behavior. 

From  the  theory  of  the  Laplace  transfonn,  we  can  get  the  original  functions  in  the 
time  domain,  after  evaluating  an  integral  of  the  fonn: 
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(48) 


p(9s ,  r)=  |  P  (s)exp(sr)ds, 

C—i  oo 


where  c  is  a  real  number  to  the  right  of  all  the  singularities  of  Ps .  For  ideal  gases,  as  is 
the  case  in  this  work,  the  right-hand  sides  of  Eqs.  (43)  and  (46)  do  not  have  real  or 
imaginary  poles  except  5  =  0.  The  other  singularities  are  the  branch  points  at  s  =  ±i ,  see 
Ref.  24.  After  closing  the  integration  contour  to  the  left,  and  defining  a  branch  cut  along 
the  imaginary  axis  between  the  branch  points  s  =  +/  and  s  =  —i ,  we  get  the  following 
integral  representation  for  the  pressure  perturbation  amplitude  at  the  shock  as  a  function 
of  time: 


Ps  (0  =  — (R  - 1)  f  cos (rsz)z2f  (. z)dz  +  [  (l  -  cos (rsz))f  ( z)dz ,  (49) 

ni  n  i 


2a , 


where  the  time  variable  rs  is  found  from  (41)  at  x  =  M2a2t 


kDtyj  l  -  M2 

RM2 


(50) 


and  the  function  /  (z)  in  (49)  is  given  by 


L(z): 


VT" 


!(l-z2)+(a10z2-an)2 


(51) 


If  we  had  imposed  different  boundary  conditions  behind  the  shock,  other  than  the  isolated 
shock  condition  studied  in  this  Section,  we  would  have  obtained  another  auxiliary 
function  /  (z)  inside  the  integral,  in  Eq.  (49).  For  example,  for  rigid  piston  boundary 

conditions,  or  free  surface  conditions  at  x  =  0 ,  the  function  /  would  be  the  solution  of 


an  inhomogeneous  functional  equation,  which  has  to  be  solved  with  techniques  similar  to 
those  outlined  in  Ref.  16,  see  Appendix. 
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The  Laplace  functions  obtained  before  are  a  very  powerful  tool.  Indeed,  it  is  not 
difficult  to  arrive  to  the  late  time  asymptotic  for  the  pressure  perturbations  at  the  shock 


front.  For  (49)  and  (51)  we  obtain: 


24 


~  ,  ,  ~oo  cos  r  -  sin  r  (  _5/2\ 
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where  the  quantity  o  s  is  defined  as 
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(53) 


and  the  constant  asymptotic  pressure  perturbation  is 


Ps  =£i 


2RM; 

y  +  1 


(54) 


Operating  in  a  similar  way  with  the  Laplace  function  of  the  shock  front  ripples 
(46),  we  get  the  exact  temporal  evolution: 

e*(Y  +  l)  f 


khxs{t)  = 


4kM2-J\-M 


2  J 

2  0 


M2(R-  1)z-^ 


sin 


irsz)f(z)dz. 


(55) 


where  rs(t)  and  /  (z)  are  defined  by  (50)  and  (51),  respectively.  At  t  —>  °°  we  find  from 


(55)  that  kbxs  approaches  the  asymptotic  value  given  by  (34),  as  it  should  be.  It  is  clear 
that  the  late  time  asymptotic  of  the  shock  ripple  would  also  show  oscillations  with  an 
amplitude  that  decays  as  fi/2 ,  superimposed  upon  the  constant  value  given  by  (34).  As 
noted  before,  the  asymptotic  amplitude  of  the  shock  front  corrugation  has  opposite  sign 
with  respect  to  that  of  the  pressure  disturbance  at  the  same  position  on  the  shock  surface 
Figure  3  shows  the  results  of  the  solution  for  a  strong  shock  with  Mx  =  10  in  an 
ideal  gas  with  y  =  5/3  interacting  with  the  “stripes”  pre-shock  density  modulation.  The 
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dimensionless  shock  ripple  amplitude  kdxs  in  Fig.  3  is  normalized  with  respect  to  the 
dimensionless  density  variation  amplitude,  £k  =  8pk  /  p , ,  so  that  its  asymptotic  value, 
according  to  (34),  is  equal  to  -  Q  .  This  value  corresponds  to  the  straight  horizontal  line. 
Its  negative  sign  means  that  asymptotically  the  shock  front  is  farther  from  the  piston 
where  the  pre-shock  density  is  initially  lower,  as  shown  in  Fig.  2.  Figure  3  shows  the 
normalized  dimensionless  shock  ripple  amplitude  vs.  the  nonnalized  time  kDt .  It 
demonstrates  that  although  the  time  history  of  the  shock  ripple  amplitude  depends  on  the 
boundary  condition  at  the  piston  surface,  its  constant  asymptotic  value  does  not,  as 
predicted  by  Eq.  (34).  This  example  also  demonstrates  that  the  shape  of  a  strong  shock 
wave  propagating  along  the  “stripes”  reaches  its  asymptotic  value  after  it  travels  slightly 
more  than  one  lateral  wavelength.  For  the  weak  shocks,  one  wavelength  is  enough.  On 
the  other  hand,  the  oscillations  of  khxs  around  its  asymptotic  value  (24)  decay  slowly,  cf. 
(52),  decreasing  to  the  -10%  level  as  the  strong  shock  propagates  -16  wavelengths.  This 
slow  decay  of  oscillations  is  probably  a  feature  of  the  ideal-gas  equation  of  state  used 
here,  cf.  Ref.  25,  and  particularly  Fig.  1  there. 

B.  Pre-shock  density  modulation  localized  in  the  lateral  direction  -  “joint”  and 
“tube” 

It  is  of  practical  interest  to  calculate  the  response  of  a  shock-piston  flow  to  a 
spatially  localized,  rather  than  periodic,  pre-shock  density  perturbation.  An  important 
particular  case  is  a  density  perturbation  concentrated  in  a  planar  layer,  perpendicular  to 
the  shock  front,  Fig.  1(d),  which  emulates  a  “joint”  in  a  laser  target  or  a  “thermal  layer” 
near  a  heated  wall.10'14  Assuming  a  cylindrical,  rather  than  planar,  symmetry  of  the 
density  perturbation  shown  in  Fig.  1(d),  we  obtain  the  case  of  a  “tube,”  a  density 
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perturbation  concentrated  in  a  cylinder  perpendicular  to  the  shock  front,  which  emulates  a 
fill  tube  in  a  laser  target.9 

To  solve  this  problem,  we  decompose  the  pre-shock  density  perturbation  into 
Fourier  integral.  For  each  of  the  Fourier  modes,  the  problem  is  reduced  to  the  “stripes” 
case  of  Fig.  1(a),  which  has  already  been  solved  above.  We  only  need  to  substitute  this 
solution  into  the  Fourier  integral  and  to  perform  the  integration. 

We  start  with  a  simpler  “tube”  case  and  limit  ourselves  to  a  Gaussian  pre-shock 
density  profile,  which  is  decomposed  into  a  Fourier  integral  as 


^  =  e  exp(-  ?  V  )  =  -  Vf  fcft  exp 

P,  47 iq  l 


(  yp  Y? 
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where  1  Iq  is  the  characteristic  diameter  of  the  pre-shock  density  perturbation  and  £  (/  «  1 
is  its  small  dimensionless  amplitude.  The  expression  under  the  integral  over  k  is  the 
Fourier  amplitude  8pA.  /  Pj  of  the  planar  periodic  density  perturbation  mode 
corresponding  to  the  wavenumber  k.  For  each  of  these  modes,  the  solution  is  given  by 


£8x^(0  =  —' <t>(kDt) 


(57) 


where  x  =  kDt  is  the  dimensionless  time  for  the  k- th  mode,  <t>(x  )  is  a  known  function, 
the  same  for  all  modes  (its  typical  graphs  are  shown  in  Fig.  3).  Since  the  shock  front  is 
initially  planar,  <t>(0)  =  0  ;  the  limit  of  this  function  at  large  argument  is  given  by  (34). 


Therefore,  at  any  time  the  shape  of  the  perturbed  shock  front  is  given  by 

1/0(A r)<t>(kDt) . 


qhxs  (r,  t )  =  —  f  dk  expl 
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For  large  t,  the  argument  of  $  in  (58)  is  large  everywhere  except  at  k  -  0 ,  and  we 
can  approximate  this  function  in  the  Fourier  integral  by  —  Q ,  its  asymptotic  value  at  large 


argument,  see  (34).  Then,  perfonning  integration  over  k,  we  find  the  time-independent 
asymptotic  shape  of  the  shock  front  for  the  “tube”  case: 


q§xs(.r,t  — >  °o)  =  —eq  -^-(9  exp  ~q2r2  h 


by  the  “tube”  kind  of  the  pre-shock  density  perturbation  always  remains  finite,  localized 
in  the  vicinity  of  the  “tube,”  at  the  distance  of  the  order  of  its  diameter,  in  qualitative 
agreement  with  the  numerical  results  of  Ref.  9. 

For  the  planar,  rather  than  cylindrical,  symmetry  of  the  localized  pre-shock 
density  perturbation,  the  “joint”  case,  the  Fourier  integral  is  written  as 
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(62) 


8p(A,~^  =  e  exp(- q 2 y 2 )  =  — j=  ? dk  exp[ - Icos ky  , 

Pi  t  4qJ 

where,  as  above,  \/q  is  the  characteristic  width  of  the  pre-shock  density  perturbation  and 
£  «  1  is  its  small  dimensionless  amplitude.  Therefore  the  shape  of  the  shock  front  is 
given  by 

e  7  f]k  (  k2  \ 

qbxs(y,t)  =  -F=\—ZK P  -—2  ®(kDt)cosky  .  (63) 

Vtt  o  k  [  4 q~  J 

The  asymptotic  estimate  for  this  “joint”  case  cannot  be  made  similarly  to  (58) 
because  Fourier  integral  (63)  diverges  at  its  lower  limit  k  -  0  if  '(kDt )  is  simply 
replaced  by  its  constant  asymptotic  value  -  Q  .  The  reason  for  this  divergence  is  simple. 
As  seen  from  (34)  and  (57),  the  asymptotic  shock  displacement  amplitude  is  large  for 
long  wavelengths  (small  k).  In  the  above  “tube”  case,  the  long- wavelength  modes  occupy 
a  small  portion  of  the  phase  space,  which  is  why  their  contribution  to  the  Fourier  integral 
is  small.  In  contrast  with  this,  for  a  “joint’  case  the  long-wavelength  modes  make  the 
dominant  contribution  to  the  Fourier  integral  (63),  making  it  diverge  at  t  — >  °o .  A  simple 
calculation  yields  the  time-dependent  asymptotic  shape  of  the  shock  front  for  the  “joint” 
case: 

q$xs  (y,t  — » °°)  =  —eq  -J= 

where  2F2(a,b;c,d;z )  is  the  generalized  hypergeometric  function, 
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and  yE  =  0.5772...  is  the  Euler’s  constant. 


(65) 
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Equations  (66)  and  (67)  predict  that  the  shock  perturbation  exhibits  two  spatial 
scales,  in  qualitative  agreement  with  the  simulation  results11  and  the  experimental 
findings.  ’  In  the  vicinity  of  the  “joint”  plane  its  shape  retains  the  characteristic  length 
scale  1  /  q  ;  in  particular,  the  radius  of  its  curvature  approaches  a  steady-state  value  of  this 
order,  which  is  independent  from  the  boundary  condition  at  the  piston,  similarly  to  the 
“tube”  case.  On  the  other  hand,  Eq.  (66)  predicts  unlimited,  albeit  slow,  logarithmic 
growth  of  the  shock  displacement  amplitude.  It  contains  the  tenn  A  that  depends  on  the 
time  history  of  the  shock  displacement  amplitude,  and  therefore,  depends  on  the 
boundary  conditions  at  the  piston.  Note  that  the  asymptotic  shape  of  the  shock  profile 
(66)  at  late  times,  qDt  »  1 ,  at  large  distances  from  the  “joint”,  |  qy  |»  1 ,  is  self-similar, 


depending  on  y  and  t  only  via  the  dimensionless  combination  Dt/\y\  and  independent  of 


q.  Figures  5(a),  (b)  demonstrate  that  this  is  indeed  the  case.  The  perturbed  shock 
dynamics  is  different  for  the  cases  of  rigid  piston  (a)  and  free  surface  (b)  because  the 
parameter  A  equals  -0.757  in  the  former  case  and  1.71  in  the  latter.  The  growth  rates, 
however,  approach  the  same  limit  for  both  cases. 
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C.  Double-periodic  pre-shock  density  modulation  -  “quilt” 

Now  consider  the  “quilt”  case  of  double-periodic  pre-shock  density/entropy 
perturbation,  Fig.  1(b),  for  the  isolated  shock  case.  By  proceeding  much  in  the  same  way 
as  in  Section  III. A,  we  get  the  following  exact,  closed  form  expressions  for  the  Laplace 
transforms.  Generalization  of  (43)  and  (46)  for  £  0  >  0  gives,  respectively, 
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The  time  evolution  of  the  pressure  and  shock  ripple  perturbations  can  also  be 
easily  obtained,  as  a  function  of  the  shock  intensity,  fluid  compressibility  and  the 
parameter  .  In  both  Eqs.  (68)  and  (69)  we  recognize  the  denominator  s2  +C,2 ,  whose 
poles  at  s  =  ±/£0  give  rise  to  oscillations  of  constant  amplitude  in  the  domain  of  the 
variable  rs,  the  frequency  of  which  is  £  0 .  This  oscillation  is  convoluted  with  another 
oscillating  function.  We  get,  for  example,  for  the  shock  front  ripple: 

kbxs ( t )  =  j[q  sin( rsz)  +  g(rs , z)]zfp (. z)dz -£k  ,  (70) 

where  the  auxiliary  function  g(rs,z )  inside  the  integral  (70)  and  the  coefficient  c,  are 
given  by: 


,  (y  +  l)«20  zsin(r^0)-^0sin(/;z) 

2nM2^-M2  Cot2  "Co) 

_e,Vl -Ml 
^  n  M\  ’ 


(71) 

(72) 


Phys.  Plasmas  14,  072706  (2007)  30 


and  /  (z)  is  defined  by  (51).  To  get  an  approximate  analytical  formula  valid  when  the 


shock  is  far  enough  from  the  surface  x  =  0 ,  we  must  distinguish  the  two  regimes:  long- 
wavelength  C0<1  and  short-wavelength  £  0  >  1 ,  as  defined  in  Section  II. 

For  the  long-wavelength  regime,  the  pressure  perturbation  at  late  time  is  given  by 
the  asymptotic  fonnula 

Ps(0  =  K  cos {kxDt)+b02  sin (kxDt),  (73) 

where  the  coefficients  b0l  and  b02  are  given  by: 
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(74) 


For  the  shock  front  ripple  amplitude  we  reproduce  the  asymptotic  formula  (32). 

Increasing  the  longitudinal  wavenumber  further,  we  enter  the  short-wavelength 
regime,  where  C,  0  >  1 .  The  asymptotic  fonnula  for  the  shock  pressure  perturbation 
amplitude  becomes: 

ps(t)  =  e0  cos(k xDt)  ,  (75) 


with 


c-o  —  I -  .  V/vV 

—  1  —  ®ioCo 

The  asymptotic  formula  for  the  shock  ripple  amplitude  is  (35). 

The  long-wavelength  condition  (26)  for  the  same  parameters  as  in  Fig.  3, 

Mj  =  10  ,  y  =  5/3,  requires  kx  / k  <  0.5049 .  Choosing  the  wavelength  ratio  in  this  range, 

kx  Ik  =  1/5 ,  we  expect  the  asymptotic  expression  (32)  to  be  valid.  Figure  6 
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demonstrates  that  this  is  indeed  the  case,  for  the  rigid  piston  and  free  surface  boundary 
conditions  as  well  as  for  the  isolated  shock.  As  in  Fig.  3,  the  time  history  of  the  shock 
ripple  amplitude  depends  on  the  boundary  condition  at  the  piston,  whereas  its  asymptotic 
behavior  does  not. 

Figure  7(a)  illustrates  the  short-wavelength  “quilt”  pre-shock  density  modulation 
with  kx/kv  =  1 ,  the  other  parameters  being  the  same  as  above,  for  isolated  shock  case. 

We  see  that  the  asymptotic  amplitude  and  phase  predicted  by  (35)  are  achieved  very 
rapidly. 

This  asymptotic  expression,  however,  is  not  applicable  in  the  presence  of  a  piston 
driving  the  shock  wave  into  the  short-wavelength  “quilt”.  In  this  case,  not  covered  by  the 
elementary  theory  of  Section  II. B,  sonic  waves  reflected  from  the  piston  can  reach  the 
shock  wave  and  affect  its  behavior.  The  appropriate  asymptotic  formulae  can  still  be 
obtained  using  the  Laplace  transform  technique.  Here  we  present  the  results,  whose 
derivation  is  outlined  in  the  Appendix. 

For  the  short  wavelength  regime  with  sufficiently  large  C,0  the  sonic  waves 
radiated  by  the  shock  at  the  nonnalized  frequency  C,  0  given  by  (25)  that  corresponds  to 
the  physical  frequency  co  =  kxD  can  be  reflected  from  the  piston  surface  and  catch  the 

shock  front  from  behind.  This  will  occur  if  the  x-velocity  of  the  reflected  sonic  wave  is 
greater  than  the  shock  speed  in  the  reference  frame  of  the  shocked  fluid.  That  is,  the  angle 
cp  that  the  sonic  ray  forms  with  the  x-axis  (nonnal  to  the  piston  surface)  must  satisfy 

coscp  >  M2 ,  see  Ref.  18.  This  condition  can  be  expressed  as 
„  ,  1  +  dfy2 

>  cosh20  = - (77) 

0  s  1  -Ml 
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Let  us  define  the  auxiliary  quantity  o0  by  =  cosho0 .  Equation  (77)  is  obviously 
equivalent  to  o0  >  20  s ,  which  is  also  seen  to  be  the  requirement  for  a  “successful 
reflection  at  the  piston  boundary”.  When  this  condition  is  met,  the  shock  starts  to  oscillate 
with  another  nonnalized  frequency  C,2  given  by: 

C2=cosh(a0-2eJ.  (78) 

The  higher-order  terms  appear  as  follows.  If  the  sonic  waves  emitted  at  the 
frequency  C,2  could  reach  the  shock  surface  from  behind,  after  being  reflected  at  the 
piston,  another  mode  (with  a  lower  frequency^)  would  be  excited  at  the  shock  front.  For 
this  to  happen,  a  requirement  similar  to  Eq.  (77)  has  to  be  satisfied,  which  is  obtained  by 
replacing  20  s  with  40  s ,  60s ,  ...,  etc.,  in  the  right-hand  side  of  Eq.  (77).  The 
corresponding  normalized  oscillation  frequencies  therefore  are 

C2b  =cosh(o  o-2/i0  J  .  (79) 


downshifted  in  the  reference  frame  of  the  shocked  fluid  (with  a  normalized  frequency 
C»,+i  <  C  )•  Then,  after  a  finite  number  of  reflections,  the  frequency  of  the  last  excited 


mode  would  be  eventually  so  low  that  the  x-velocity  of  the  associated  sonic  ray  is  less 
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than  the  shock  speed  with  respect  to  the  shocked  fluid,  a2M2 ,  and  then  the  process  of 
“successful  reflection”  terminates.  This  phenomenon  is  totally  equivalent  to  the  problem 
of  spontaneous  acoustic  emission  (SAE)  for  a  shock  moving  into  a  certain  non-ideal  gas, 
in  the  presence  of  a  reflecting  boundary  from  behind,  as  discussed  in  Ref.  24.  The 
difference  with  the  problem  being  studied  here  is  that  in  our  case,  the  shock  is  forced  to 
oscillate  because  of  the  pre-shock  density  modulation,  so  that  the  oscillations  take  place 
for  arbitrary  EOS,  including  ideal  gas.  In  Ref.  24,  on  the  contrary,  the  initial  oscillation 
was  provided  by  the  SAE  phenomenon,  due  to  the  non-ideal  EOS  of  the  fluid  into  which 
the  shock  was  moving.  The  interesting  point  to  note  is  that,  once  we  make  the  shock  to 
oscillate  with  any  initial  frequency  (asymptotically  in  time,  independently  of  the  reason 
of  that  oscillation),  all  the  other  additional  modes  that  are  also  seen  to  be  excited  at  its 
surface,  are  just  due  to  this  process  of  successful  reflection  downstream,  and  their  origin 
is  therefore  essentially  kinematic. 

The  asymptotic  pressure  perturbation  and  shock  front  ripple  amplitudes  are 
expressed  as 


i  v 

Ps(t)  =  Y<e2nRnp  C0S 


n=  0 


r  3 
y^KDt 

S  o 


(81) 


khxs(t)  =Y^j2nRnp  sin 


n= 0 


r  7 
y^KDt 

S  o 


(82) 


Here  Rp ,  as  above,  is  the  reflection  coefficient  for  the  sonic  waves  incident  upon  the 
piston.  As  in  (38)  and  (40),  this  factor  accounts  for  the  phase  and  amplitude  change  of  the 
sonic  wave  at  reflection  from  the  piston.  The  coefficient  e0  is  given  by  (76);  the  other 


coefficients  are 
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The  value  of  /„  given  by  (84)  is,  of  course,  the  same  as  the  coefficient  at  sin (kxDt)  in 
(35).  For  the  isolated  shock  case  Rp  — »  0  and  (81)  and  (82)  are  reduced  to  (75)  and  (35), 

respectively.  Otherwise  the  terms  in  the  sums  (81),  (82)  decrease  with  increased  n  very 
rapidly  because  the  reflection  coefficient  of  the  sonic  waves  from  the  shock  front  is  lower 
for  oblique  incidence  than  for  nonnal  incidence  (39),  and  for  each  higher-order  wave 
returning  to  the  shock  front  from  the  piston  the  angle  of  incidence  is  more  oblique  than 
for  the  previous  one.  Therefore  in  most  cases  the  two  first  terms  in  the  sums  (81),  (82) 
provide  sufficient  accuracy.  In  the  limit  kx  /k  — »  °° ,  C,  0  —>  °°  we  have 


C2«/^o  ~ ^  [(1-M2)/(1  +  M2)]\  and  (82)  is  reduced  to  (40),  as  it  should  be. 

When  a  piston  driving  the  shock  wave  is  present,  the  lower- frequency  oscillations 
superimposed  upon  the  oscillations  at  the  driving  frequency  co  =  kxD  are  clearly  seen  in 
Fig.  7(b).  The  asymptotic  fonnula  (82)  is  shown  to  provide  good  accuracy  with  just  one 
extra  term,  corresponding  to  n  —  1 ,  retained. 
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D.  Pre-shock  standing  sonic  wave 

All  the  above  results  refer  to  the  density  modulation  ahead  of  the  shock,  which,  at 
constant  pre-shock  pressure,  corresponds  to  an  “entropy  wave”  that  is  stationary  with 
respect  to  the  resting  fluid  ahead  of  the  incident  shock  front.  In  general,  a  small 
perturbation  of  a  unifonn  resting  fluid  is  decomposed  into  a  density/entropy  perturbation, 

1  o 

a  vorticity  perturbation  and  a  sound  wave.  Our  theory  can  be  easily  modified  for  the 
case  when  the  perturbation  ahead  of  the  shock  front  is  either  a  standing  sonic  wave  or  a 
periodic  vortex  system.  We  consider  here  the  fonner  case  as  an  example. 

We  limit  ourselves  to  the  situation  when  the  wave  vectors  in  the  two 
superimposed  running  sonic  waves  that  make  a  standing  wave  are  parallel  to  the  shock 
front,  in  the  y  direction.  (This  is  because  there  is  no  reflection  of  the  incident  sonic  waves 
from  a  shock  front,  so  we  cannot  form  any  other  standing  waves.)  Using  the  above  tenns, 
the  perturbation  pattern  at  any  instant  looks  like  “stripes”,  Fig.  1(a).  The  areas  of  high  and 
low  density,  however,  change  places  periodically,  which  is  why  the  shock  front  feels  this 
upstream  perturbation  as  a  “quilt”,  Fig.  1(b),  with  pre-shock  density  and  other  variables 
oscillating  at  y  =  0  .  For  a  given  lateral  wavenumber  k ,  the  frequency  of  these  pre-shock 

oscillations  is  0),  =  kai ,  where  ax  is  the  pre-shock  speed  of  sound,  see  Fig.  2.  The 
perturbations  coming  to  the  shock  front  at  the  time  t  are 

5P*(0  =  e*PiCos(ov),  dpk(t)=£kplcos{alt), 

(  \  (85) 

dvxk( t)  =  0,  dvyk( t)  =  -/£,«,  sin(oy) , 

where  the  dimensionless  amplitude  of  the  sonic  wave  is  £k  «  1 .  The  density 
perturbation  (85)  is  similar  to  that  in  a  “quilt,”  Fig.  1(b),  with  the  effective  longitudinal 
wavenumber 
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(86) 


The  perturbations  of  the  pressure  and  lateral  velocity  are  not  present  in  the  “quilt.”  To 
account  for  these,  the  boundary  conditions  ( 16)-(  19)  need  to  be  modified  appropriately. 

Obviously,  for  a  sufficiently  strong  shock  wave,  the  ratio  kx  eff  /  k  becomes  small 

enough  to  satisfy  (26),  and  we  approach  the  long-wavelength  “quilt”  case.  Physically,  the 
post-shock  speed  of  sound  is  proportional  to  the  shock  velocity,  and  the  pre-shock 
oscillation  frequency  C0j  =  ka]  for  sufficiently  strong  shocks  is  much  smaller  than  the 


corresponding  post-shock  frequency  to  =  ka2 .  The  resulting  requirement  on  the  shock 
strength  is 


Ml  > 


(87) 


When  the  shock  wave  is  sufficiently  strong  the  similarity  between  the  cases  of  the 
upstream  standing  sonic  wave  and  the  “quilt”  increases  approaching  equivalence.  This  is 
because  both  velocity  and  pressure  behind  the  strong  shock  wave  are  so  much  higher  than 
their  pre-shock  values,  that  the  pre-shock  perturbations  of  these  variables  can  be 
neglected.  On  the  other  hand,  density  variation  in  the  shock  wave  is  limited,  hence  the 
perturbed  strong  shock  wave  mostly  reacts  to  the  pre-shock  density  perturbations. 

Figure  8  shows  the  time  histories  of  shock  ripple  amplitude  for  y  =  5/3  and  two 

Mach  numbers,  M,  =  10  and  25.  For  y  =  5/3  Eq.  (87)  gives  Mx  >  V3  =  1.732  ,  hence  both 
of  these  cases  satisfy  the  long-wavelength  condition.  The  agreement  with  (32)  is  good,  as 
it  should  be. 
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E.  Richtmyer-Meshkov-like  instability  and  lateral  mass  redistribution 
1 .  RM-type  perturbation  growth  at  the  inner  interface  and  free  surface 

It  has  been  demonstrated  in  Ref.  1  that  a  fast  application  of  an  external  pressure  to 
a  rippled  surface  produces  a  fingering  instability  similar  to  the  classical  RM  instability, 
RM-type  instability.  We  show  here  that  such  instability  develops  even  if  the  loaded 
surface  or  interface  is  initially  planar.  If  the  half-space  behind  the  surface  x  =  0  contains 
pre-shock  non-uniformities,  then  the  perturbations  needed  to  drive  the  instability  are 
instantly  deposited  at  this  surface,  when  a  planar  shock  wave  either  reaches  this  surface 
(isolated  shock  case)  or  is  launched  from  it  (all  other  cases  involving  a  non-rigid  piston) 
and  brought  to  the  surface  by  sonic  waves. 

After  t  =  0+  ,  the  free  surface  x  =  0  becomes  a  vortex  sheet,  which  exhibits  a 
linear  ripple  growth  as  in  the  Richtmyer-Meshkov  instability. 1  In  the  case  of  isolated 
shock  wave,  this  surface  becomes  a  small-amplitude  weak/contact  discontinuity 
separating  two  regions  of  fluid:  to  the  left,  x  <  0 ,  we  would  have  irrotational  velocity 
field  generated  by  the  left-facing  sonic  waves  radiated  by  the  rippled  shock  front  since 
t  =  0+  .  To  the  right,  at  0  <  x  <  M2a2t  in  the  frame  of  reference  of  the  shocked  mass,  the 
fluid  is  filled  with  sonic  waves  plus  vorticity  and  entropy  perturbations  generated  at  the 
shock  front.  We  show  here  that  this  interface  exhibits  a  RM-type  linear  growth,  too. 

The  Laplace  transform  technique  makes  it  possible  to  find  the  asymptotic 
behavior  of  the  x  =  0  surface  for  these  cases.  For  the  “stripes”  and  long-wavelength 
“quilt”  cases  of  pre-shock  density  modulation,  the  velocity  perturbation  amplitude  at  the 
surface  x  =  0  tends  to  a  constant  asymptotic  value 


v  = 


-0+ 


[coshe^  +  a^ejjRxej 
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20 


sinh  (0s)+Co 
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where  the  velocity  perturbation  at  t  =  0+  is  given  by  vxi+=-M2(R- 1)/2 ,  see  (17).  For 
the  isolated  shock  case  the  function  Ps  (0  )  defined  by  (42)  is  explicitly  determined  by 
(43)  or  (68).  The  procedure  of  calculating  its  value  for  the  free  surface  case  is  outlined  in 
the  Appendix,  where  the  function  OCj(0s)  is  defined  by  Eq.  (A3). 

For  the  short- wavelength  “quilt”,  sonic  waves  radiated  downstream  by  the  rippled 
shock  front  reach  the  x  =  0  surface  and  make  it  oscillate.  The  number  of  the  oscillating 
modes  is 


N'  ■ 


In  f.  +  VfH)  1 

i„E^  2 

1  -M, 


(89) 


cf.  (80).  This  number  equals  either  N  given  by  (80)  or  N  + 1 ,  because  a  sonic  wave 
radiated  by  the  shock  front  that  reaches  the  piston  is  not  necessarily  able  to  catch  up  the 
shock  front  after  its  reflection.  The  dimensionless  frequencies  of  these  modes  are 

C 2n-\  =  cosh(o 0  (2/?  1)0 5 ) ,  /I  =  1,2,  ...,  N' ,  (90) 

cf.  (79).  For  the  isolated  shock  case,  only  the  n  =  1  mode  can  be  excited  at  the  x  =  0 
surface  because  there  is  no  reflection  of  sound  from  this  surface  back  to  the  shock. 


The  asymptotic  expression  for  the  velocity  perturbation  is 


N ' 

v^v~+lviM_x  cos 


C  2h-i  kxDt 


n—\ 


Co  -yjl  —M 


(91) 


where  v”  is  given  by  (88)  and  the  oscillation  amplitudes  are 
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For  the  isolated  shock  case,  TV'  in  (91)  is  replaced  by  the  smaller  of  AT'  and  1. 
We  see  that  the  ripples  at  the  free  surface  grow  because  of  the  RM-type  instability 
developing  on  that  surface.  Superimposed  on  this  linear  growth  are  the  oscillations 
induced  by  the  waves  radiated  by  the  shock  surface,  which  is  forced  to  beat  with  a 
rhythm  imposed  by  the  longitudinal  density  modulation  ahead  of  it.  Similar  arguments 
would  apply  to  the  case  of  a  rigid  boundary  at  x  =  0  .  The  difference  is  that  for  a  rigid 
piston  it  is  the  sonic  pressure  and  tangential  velocity  perturbations  at  x  =  0  which  would 
show  the  oscillations  with  the  nonnalized  frequency^, .  In  the  limit  K  /  ky  OO  9 

C,0  — >  °o,  C 2,-1 /Co  ~ >  [(1  -M2)/(l  +  M2)]"~1'2 ,  and  the  sum  over  n  in  (91)  is  reduced  to 
the  solution  (37)  and  (38). 

Figure  9  is  plotted  for  an  isolated  shock  wave  interacting  with  a  “quilt”  pre-shock 
density  modulation,  Fig.  1(b).  In  the  long-wavelength  case,  when  (26)  is  satisfied,  the 
periodic  shock  oscillations  do  not  generate  running  waves  that  can  reach  the  x  =  0 
surface.  Therefore  the  longitudinal  velocity  perturbation  amplitude  at  x  =  0  approaches 
the  constant  asymptotic  value  given  by  Eq.  (88),  similarly  to  all  other  examples  of  the 
RM-type  instability.1' 16  It  is  illustrated  by  Fig.  9(a)  plotted  for  kjky  =  1/5 ,  y  =5/3, 

Mj  =  10  ;  the  x-velocity  perturbation  amplitude  V  is  normalized  with  respect  to  £,  D  . 
The  constant  asymptotic  value  of  V ’  means  that  the  initially  planar  contact  interface 
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formed  by  the  fluid  particles  initially  located  at  x  =  0  becomes  rippled,  with  the  ripple 
amplitude  growing  linearly  with  time,  which  is  typical  for  the  RM-type  instabilities.1 
Figure  9(b)  demonstrates  a  similar  calculation  for  a  short-wavelength  “quilt,”  kx  / k  =  1 . 

In  this  case,  the  oscillating  shock  front  radiates  sonic  waves  that  reach  the  x  =  0  surface. 
As  a  result,  its  velocity  V  oscillates  around  the  constant  average  value  given  by  the  same 

Eq.  (88),  rapidly  approaching  the  asymptotic  harmonic  oscillation  predicted  by  Eq.  (91). 
In  this  case,  too,  the  ripple  amplitudes  of  the  contact  interface  grow  linearly  with  time,  on 
the  average,  with  a  non-decaying  oscillation  superimposed  upon  the  linear  growth.  This 
growth,  however,  is  not  easy  to  detect  in  either  case,  even  in  a  simulation,  unless  it  is 
done  with  a  Lagrangian  code  that  tracks  individual  fluid  particles.  Indeed,  the 
unperturbed  density  is  equal  on  both  sides  of  the  rippled  interface,  and  therefore  its 
growth  is  not  captured  by  any  diagnostic  that  tracks  lateral  mass  redistribution,  whether  it 
is  numerical  or  experimental,  like  face-on  radiography. 

2.  Lateral  mass  redistribution 

The  lateral  mass  redistribution  is  a  very  important  manifestation  of  all  the 
interfacial  instabilities,  including  RT,  classical  RM,  ablative  RM  and  other  RM-type 
instabilities,  which  allows  us  to  directly  observe  their  development  in  laser  targets  with 
the  aid  of  face-on  x-ray  radiography,  as  in  Refs.  26-28.  Our  next  example  of  the  RM-type 
growth  that  involves  the  lateral  mass  redistribution  refers  to  the  “stripes”  case  of  pre¬ 
shock  density  modulation,  Fig.  1(a),  with  the  free  surface  (constant  pressure)  boundary 
condition  at  the  piston.  The  free  surface  is  initially  planar.  However,  the  ripples  rapidly 
develop  at  the  surface  because  the  shock  wave  driven  by  a  constant  pressure  propagates 
faster  into  the  lighter  fluid,  and  the  piston  follows  it.  Then  these  ripples  continue  their 
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linear  growth,  similarly  to  all  the  RM-type  situations.1  Figure  10(a)  shows  the  time 
history  of  the  nonnalized  ripple  amplitude,  kbxp  and  of  its  velocity  normalized  as  in  Fig. 

9,  V  =d[khx p  )/d(kDt ) .  We  see  a  good  example  of  a  RM-type  behavior  -  linear 
asymptotic  growth,  constant  velocity  given  by  (88).  Just  as  above,  the  growth  is  in  the 
negative  direction  because  the  piston  lags  behind  where  the  density  is  higher,  Sp^.  >  0  ,  see 
Fig.  2. 

For  our  piston  problem,  the  areal  mass  modulation  amplitude  8m  can  be 
expressed  as 


8m(t)  = 


p2X 


2% 


■  k8x  +  — — -  k8x  +k\ 

p  R  l 


8p(x,t)  (5P](0 


v  p2  Pl  j 
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(93) 


where  the  time-averaged  value  of  the  pre-shock  density  modulation  at  the  shock  front  is 
defined  as 


(Sp  j  (0)  =  -  JSp  j  (x  =  Dt')dt' ;  (94) 

e.  g.,  for  the  “stripes”  case  (Sp, )  /  p,  =ek  =  const ,  whereas  for  the  “quilt”  case  (20) 

(Spj)/  Pj  =ek  sin (kxDt)/(kxDt),  cf.  Refs.  1,  27,  28.  The  three  contributions  to  5 m  in  the 

square  brackets  come  from  the  ripples  at  the  piston  surface  and  the  shock  front  and  from 
the  density  modulation  left  between  the  shock  and  the  piston,  respectively.  The 
expression  (93)  accounts  only  for  the  time-dependent  areal  mass  modulation  resulting 
from  the  lateral  mass  redistribution,  it  does  not  include  the  constant  initial  mass 
modulation  amplitude,  which  for  the  “stripes”  case  of  Fig.  1(a)  equals  SpA.Z  ,  where  L  is 

the  target  thickness.  According  to  (93),  the  negative  growth  of  Sx/;  shown  in  Fig.  9(a) 
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adds  a  positive  contribution  to  dm  ,  as  it  should  be  -  the  bubble  deepens  where  the 
density  was  initially  low,  thereby  increasing  the  initial  density  contrast.  In  all  situations 
of  RM  and  RM-type  growth,  this  contribution  to  dm  from  the  unstable  surface/interface 
rapidly  becomes  dominant,  cf.  Figures  1  of  Refs.  27  and  28.  Figure  10(b),  however, 
shows  that  this  is  not  the  case  here:  at  late  time,  we  see  a  constant,  finite  difference 
between  (d  /  dt)dm  (dm  is  normalized  with  respect  to  p2AV  2n  ,  kDt  is  the  normalized 
time)  and  the  constant  asymptotic  contribution  from  the  RM-type  unstable  piston  surface, 
-  V .  This  difference  cannot  come  from  the  shock  front  contribution,  because  the  shock 

front  displacement  amplitude  oscillates  around  a  constant  asymptotic  value,  as  shown  in 
Fig.  3,  hence  its  time  derivative  is  zero  on  the  average.  Therefore  it  is  due  to  a  lateral 
mass  redistribution  in  the  volume  of  the  shocked  target.  Since  the  difference  is  negative, 
the  perturbed  lateral  flow  moves  the  mass  from  where  the  density  was  initially  higher  to 
where  it  was  initially  lower. 

To  elucidate  the  cause  of  this  redistribution,  consider  the  rigid  piston  boundary 
condition  that  eliminates  any  contribution  of  the  perturbed  surface  to  the  areal  mass 
modulation  amplitude  dm  .  Figure  1 1(a)  shows  the  time  history  of  dm  for  this  case.  The 
longer  the  shock  wave  propagates  through  the  “stripes”  density  perturbation  field,  the 
more  the  initial  amplitude  8pA,Z  is  reduced  by  the  post-shock  flow.  Figure  1 1(b)  displays 

the  profiles  of  relative  density  perturbation  amplitude,  p  =  Sp  /  p2  ,  nonnalized  with 
respect  to  ek,  at  five  instants  of  time.  On  the  average,  p  remains  below  unity  (its  pre¬ 
shock  value  shown  by  the  horizontal  dotted  line),  so  the  corresponding  contribution  to  the 
right-hand  side  of  (93)  is  negative.  One  can  say  that  a  shock  wave  helps  the  shocked  fluid 
to  reach  greater  unifonnity  by  reducing  the  pre-shock  density  contrast  in  it;  the  shock 
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brings  more  pressure  where  the  density  is  higher,  and  these  areas  expand.  Note,  however, 
that  in  the  immediate  vicinity  of  the  shock  front  p  exceeds  unity,  and  the  direction  of  the 
lateral  mass  flow  is  different  there.  Indeed,  the  shock  front  is  convex  where  the  density  is 
initially  lower  and  concave  where  it  is  initially  higher.  It  implies  that  the  mass  near  the 
shock  front  flows  laterally  to  where  the  density  was  initially  higher,  thereby  adding  a 
positive  contribution  to  8m  .  The  flow  in  the  opposite  direction  dominates  in  the  vicinity 
of  the  piston,  and  on  balance  the  contribution  reducing  the  post-shock  density  non¬ 
uniformity  prevails. 

What  if  the  shock  wave  is  driven  by  a  laser,  with  an  ablation  front  serving  as  a 
piston?  In  the  case  of  a  uniform  target  with  front  surface  ripples,  its  instant  shock  loading 
produces  the  so-called  ablative  RM  instability:  Instead  of  a  linear  growth  similar  to  that 
illustrated  by  Fig.  10(a),  both  the  ablation  front  displacement  and  areal  mass  modulation 
amplitudes  exhibits  slowly  decaying  oscillations.  Such  oscillations  of  8m  have  recently 
been  observed  experimentally.-  ’  ’  Linear  RM-type  perturbation  growth  at  the  ablation 
front  is  suppressed  by  the  “rocket  effect”19  that  is  incorporated  in  the  boundary  conditions 
(13),  (14)  -  for  details,  see  Refs.  1,  27,  29.  For  our  case  of  the  pre-shock  “stripes”  density 
perturbation,  Fig.  1(a),  the  time  histories  of  the  nonnalized  displacement  amplitude  of  the 
ablation  front,  k8xa ,  and  of  the  areal  mass  modulation  amplitude,  8m  ,  are  shown  in 

Fig.  12  plotted  for  the  parameters  of  the  ablation  front  rD  =  0.0625 ,  M a  =  0.05  ,  which 
roughly  correspond  to  a  planar  solid  DT  target  with  X  =  24  pm  ripples  on  its  front 
surface  irradiated  with  a  0.35  pm  laser  radiation  at  3x10  W/cm  ,  see  Ref.  3 1 . 

We  see  in  Fig.  12(a)  that  the  slowly  decaying  oscillations  of  the  ablation  front 
occur  around  a  nonzero  average.  This  happens  because  at  early  time  the  ablation  front, 
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following  the  shock  front,  makes  a  dent  in  the  target  where  its  density  was  initially  lower. 
Later,  the  dent  remains,  since  the  plasma  in  its  vicinity  stays  non-uniform.  The 
contribution  of  the  rippled  ablation  front  to  the  mass  modulation  amplitude  8m  does  not 
decay,  nor  is  dominant,  as  in  the  case  of  a  rippled  uniform  target.  ’  ’  ’  Rather,  Fig. 
12(b)  demonstrates  that  the  mass  flow  in  the  volume  reducing  the  density  contrast, 
dominates  here,  similarly  to  the  case  of  a  rigid  piston,  cf.  Fig.  1 1(a). 

3.  Pre-shock  density  modulation  localized  near  the  shocked  surface 

One  of  the  issues  that  could  be  addressed  studying  shock  propagation  in  a  fluid 
with  weakly  modulated  density  refers  to  numerical  simulations  of  the  RM-type  unstable 
flows.  In  laser  fusion  studies,  one  often  has  to  model  a  very  small  roughness  of  the 
irradiated  surface.  It  is  not  easy  to  resolve  very  small-amplitude  ripples  on  a  numerical 
grid.  A  typical  remedy  for  it  is  to  consider  a  planar  surface  and  introduce  a  density 
perturbation  in  the  volume.  For  example,  a  sine -wave  ripple  of  amplitude  5x0  and  lateral 

wavenumber  k  on  the  front  surface  of  the  target  whose  density  is  p,  can  be 
approximated  by  a  density  profde  at  x  >  0 : 

P  (*> Zl  =  1  +  §x  / (x)  cos(Ay) ,  (95) 

Pi 

where  / (x)  is  a  dimensional  function  ([/]  =cnf1)  localized  near  the  surface,  which 
satisfies  the  nonnalization  condition 

J  /  (x)dx  =  1 .  (96) 

o 
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This  condition  ensures  that  the  initial  areal  mass  modulation  amplitude  in  this  flat  target 
with  inner  density  modulation  equals  p,&c0 ,  the  same  as  in  a  uniform  target  of  density  p, 

with  a  single-mode  surface  ripple  whose  amplitude  is  5x0 . 

For  example,  we  can  consider  profiles  peaked  at  the  surface  ( n  —  0  )  or  near  it 
( n  >  0 )  and  exponentially  decaying  far  from  it 

fix)  =  —iqx)"  exp  (-qx) ,  n>  0 ,  (97) 

n\ 

or  Bessel-function  profiles,  exhibiting  a  slow  oscillatory  decay  with  increased  x: 

f(x)  =  qJ0(qx).  (98) 

All  the  profiles  (97)  and  (98)  satisfy  the  normalization  condition  (96).  Increasing 
the  parameter  q,  one  can  make  these  perturbations  more  and  more  localized  near  the  front 
surface,  which  is  compatible  with  relative  smallness  of  density  perturbation  provided  that 
the  equivalent  ripple  amplitude  5x0  is  sufficiently  low.  Of  course,  the  problem  is 
detennined  by  the  dimensionless  ratio  q/k  .  Apparently,  in  the  limit  q / k  — >  °° ,  the 
density  perturbation  becomes  localized  near  the  surface  on  the  length  scale  much  less 
than  X  -  2k  /  k  ,  the  wavelength  of  the  transverse  ripple.  The  questions  are  whether  in  this 
limit  the  internal  density  modulation  is  equivalent  to  the  surface  ripple  as  far  as  the  later¬ 
time  perturbation  dynamics  is  concerned,  and  if  it  is  then  how  large  should  the  ratio  q/k 
be  so  that  we  can  make  use  of  this  equivalence  in  the  simulations. 

Instantly  loading  a  free  rippled  surface  with  a  constant  pressure,  we  generate  a 
RM-type  instability  at  this  surface.  A  similar  effect  arises  from  loading  a  planar  surface 
with  a  single-mode  density  perturbation  near  the  loaded  surface.  Let  us  compare  the 
results.  Figure  13  shows  the  nonnalized  growth  rate  of  the  areal  mass  variation  amplitude 
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(in  units  of  p,Sx0)  and  its  time  derivative  vs.  normalized  time,  kDt ,  for  the  n  =  0 ,  purely 
exponential  (a),  and  Bessel-function  (b)  profiles  of  the  areal  density  perturbation  near  the 
surface.  For  both  cases,  the  ratio  ql  k  varies  in  the  range  from  1  to  4,  y  =  5/3 ,  Mj=10. 

We  see  that  the  convergence  to  the  surface  ripple  case  with  increasing  q Ik  is 
faster  for  the  Bessel-function  profile  compared  to  the  exponential.  This  is  because  the 
former  concentrates  a  larger  part  of  the  (positive)  density  perturbation  in  the  vicinity  of 
the  surface  compared  to  the  latter.  At  larger  distance  the  density  perturbation  becomes 
negative,  but  it  turns  out  to  have  less  effect  on  the  surface  RM-type  growth. 

IV.  NUMERICAL  SIMULATIONS 

The  theoretical  results  presented  above  are  exact  and  therefore  do  not  require 

4  31  32 

numerical  verification.  Rather,  as  in  our  earlier  work,  ’  ’  "  they  can  be  helpful  in  testing 
the  accuracy  of  the  codes. 

Our  simulations  were  performed  in  two  dimensions  using  the  Eulerian  fast2d 
hydrocode  developed  at  the  Naval  Research  Laboratory  (more  details  and  further 
references  are  given  in  Ref.  31).  Two  profiles  of  laterally  periodic  pre-shock  density 
modulation  were  simulated:  the  “stripes”  case  of  Fig.  1(a),  and  the  density  modulation 
(97)  localized  near  the  target  surface.  All  our  simulations  were  done  for  a  strong  shock 
wave  Mj  =  10  ,  driven  into  an  ideal  gas  y  =  5/3  by  a  rigid  piston,  whose  velocity  was 
V  =  4x  106  cm/s.  The  lateral  perturbation  wavelength  in  all  cases  was  taken  X  =  30  pm. 

Dimensions  of  the  simulation  box  were  one  full  wavelength  X  =  30  pm  in  the  lateral 
( v)  direction  and  200  pm  in  the  longitudinal  (x)  direction.  The  rectangular  grid  was 
fixed  in  space  and  time;  its  cell  size  was  0.5  pm  in  the  longitudinal  direction  and  0.9375 
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jam  in  the  lateral  direction,  which  corresponds  to  32  zones  per  wavelength  (about  as  much 
or  more  than  one  can  afford  in  multi-mode  inertial  confinement  fusion  (ICF)  relevant 
simulations34). 

Figure  14  compares  numerical  simulation  results  and  the  theory  for  the  “stripes” 
case.  The  relative  density  modulation  amplitude  £  =  Sp  /  p,  is  varied  between  5x10  4 

and  5x  10  2 .  The  areal  mass  modulation  amplitude  is  normalized  with  respect  to  its  initial 
value,  S/;z0  =  ZSp  ,  where  L  =  200  pin.  The  theoretical  line  in  Fig.  14  is  basically  the 

same  as  in  Fig.  1 1(a),  differing  only  in  nonnalization  of  8m  (vertical  axis)  and  in  time 
units  (horizontal  axis).  Figure  14(a)  presents  the  results  obtained  with  a  time  step  initially 
chosen  A t0  =0.1  ps.  This  appears  to  be  a  reasonable  choice,  given  that 

V  A t0  /  X,  =  1 .3  x  1 CT4 .  In  the  simulations  of  similar  problems  done  for  the  ICF-relevant 

conditions,  ’  ’  where  the  laser  intensity  increases  gradually,  the  code  automatically 
adjusts  its  time  step  to  fit  the  problem.  In  our  present  problem,  the  shock  wave  of  full 
strength  is  generated  instantly.  The  code  was  not  originally  written  for  introducing  instant 
shocks,  so  that  the  initial  step  needs  to  be  prescribed.  We  see,  however,  that  our 
simulation  results  deviate  from  the  exact  theoretical  results  from  the  very  beginning.  The 
deviation  is  not  healed  with  the  passage  of  time,  as  the  shock  wave  propagates  in  the 
longitudinal  direction  for  up  to  7  lateral  wavelengths.  Moreover,  the  deviation  between 
the  exact  solution  and  the  numerical  results  grows  linearly  with  time.  Figure  14(a) 
illustrates  that  this  deviation  cannot  be  attributed  to  nonlinear  effects  due  to  finite 
amplitude  of  the  pre-shock  density  modulation:  reduction  of  the  parameter  8  by  two 
orders  of  magnitude  does  not  bring  the  numerical  result  closer  to  the  theory. 
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What  matters  here  is  a  proper  choice  of  the  initial  time  step.  Taking  the  initial 
time  step  A t0  =  1  fs  and  allowing  it  to  increase  no  more  than  by  5%  per  time  step,  we 
obtain  the  results  shown  in  Fig.  14(b).  This  choice  is  seen  to  make  the  agreement  with  the 
theory  much  better;  in  particular,  the  average  growth  rate  of  8m  is  modeled  very 
accurately  for  both  values  of  the  initial  density  modulation  amplitude,  2.5x10  3  and 
10  2 ,  even  though  the  time  step  after  the  first  1  ns  becomes  the  same  as  in  the  simulations 
illustrated  by  Fig.  14(a). 

Figure  15  shows  a  similar  comparison  for  the  localized  density  modulation  (95) 
and  (97).  Here,  we  have  chosen  q  =  k  =  2n  /(30  pm)  and  did  the  simulations  for  the 
parameter  n  in  (97)  varied  from  0  to  4.  For  all  the  profiles,  the  same  value  of  8x0  =0.1 
pm  in  (95)  was  chosen,  corresponding  to  the  initial  value  of  mass  modulation  amplitude 
8 m0  =  p  ,8x0 ;  areal  mass  modulation  amplitudes  are  shown  nonnalized  with  respect  to 
this  value.  The  maximum  value  of  the  relative  density  modulation  amplitude  in  the 
profile  (97)  equals  (2k8x0 /X)xn"e~n  / n\ ,  and  for  n  increased  from  0  to  4  it  decreases 

from  2x10  2  to  4xl0~3,  which,  as  confirmed  by  Fig.  14(b),  is  sufficiently  small  for  the 
linear  theory  to  be  applicable.  Our  theory,  as  demonstrated  by  Fig.  15(a),  predicts  for  all 
cases  oscillations  around  an  average  value,  which  is  less  than  unity.  This  is  consistent 
with  our  previous  results  (Fig.  1 1):  the  shock  passage  decreases  the  average  lateral  non¬ 
uniformity  of  the  density  distribution. 

Here  again,  Fig.  15(b)  shows  that  the  simulation  started  with  a  too  large  time  step 
A t0  =0.1  ps  yields  a  result,  which  is  qualitatively  incorrect,  and  its  deviation  from  the 
exact  solution  grows  linearly  with  time.  The  choice  of  a  very  small  initial  time  step 
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A t0  =  1  fs,  helped  to  improve  the  agreement,  see  Fig.  15(c).  However,  even  in  the  latter 

case  our  numerical  modeling  displays  a  slow  negative  linear  growth  of  dm,  which  should 
not  be  there:  it  is  a  numerical  artifact. 

V.  CONCLUSIONS 

We  have  developed  a  small-amplitude  theory  describing  the  interaction  of  planar, 
constant-strength  shock  waves  with  pre-shock  non-uniformities  distributed  in  the  volume 
of  a  shocked  fluid.  This  theory  describes  the  distortion  of  the  shock  front  and  the  post¬ 
shock  flow,  including  the  lateral  mass  redistribution.  Predictions  of  our  theory  are  in 
good  qualitative  agreement  with  the  known  results  of  numerical  simulations.  For 
example,  a  shock  wave  was  shown  to  “mix”  low-  and  high-density  mass,  thus  reducing 
the  average  density  non-uniformity  in  a  shocked  volume,  in  an  agreement  with  the 
simulations  of  Ref.  7,  where  a  foam-like  medium  was  modeled.  For  a  “joint”  pre-shock 
density  perturbation,  localized  in  one  transverse  direction,  the  leading  edge  of  the 
distorted  shock  wave  was  found  to  remain  localized  near  the  “joint”  and  maintain  a 
constant  curvature/length  scale  determined  by  the  joint  width,  whereas  the  shock  front  as 
a  whole  approaches  a  self-similar  shape,  in  agreement  with  the  numerical  results  of  Ref. 
8,  where  a  laser  target  containing  an  actual  joint  was  modeled.  A  slightly  simpler 
behavior  was  found  for  a  “tube”  pre-shock  density  perturbation,  localized  in  both 
transverse  directions:  the  shock  front  distortion  is  finite  and  contained  in  the  vicinity  of 
the  “tube,”  whereas  the  RM-type  growth  of  the  areal  mass  modulation,  not  being  limited, 
remains  localized,  in  agreement  with  the  numerical  results  obtained  in  Ref.  9  for  laser 
targets  with  feeding  tubes. 
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The  importance  of  our  results  stems  from  the  fact  that  these  details  of  shock 
distorting  and  post-shock  behavior  have  been  quantified  and  described  by  exact  analytical 
solutions.  Such  solutions  are  very  helpful  for  testing  the  accuracy  of  the  hydrodynamic 
simulation  codes,  particularly  when  such  codes  must  be  used  to  model  the  growth  of 
perturbations  starting  from  their  very  small  initial  amplitudes.  Small-amplitude 
perturbations  are  notoriously  hard  to  resolve,  and  there  is  always  a  risk  either  to 
contaminate  the  simulation  results  with  a  numerical  noise  or  to  damp  them  heavily  with  a 
numerical  dissipation/  Our  numerical  examples  with  small  pre-shock  density 
perturbations  demonstrate  that  it  does  not  take  much  for  the  numerical  solution  to  deviate 
from  the  exact  analytical  solution,  and  that  in  the  absence  of  acceleration  this  numerical 
error  tends  to  grow  linearly  with  time,  in  a  RM-type  fashion.  Unless  the  code  has  been 
thoroughly  tested  in  this  regime,  such  deviations  are  likely  to  remain  undetected  in 
realistic  multi-mode  simulations,  where  the  exact  solution  is  unknown. 

For  the  constant-strength  shocks  propagating  into  a  uniform  (on  the  average) 
density  target,  the  theory  is  essentially  complete.  It  is  very  important  to  generalize  it  for 
shock  waves  whose  strength  decreases  with  time,  including  those  propagating  up  the 
density  gradients  in  the  laser  targets.  This  advancement  is  needed  make  the  theory 
directly  applicable  to  the  situations  relevant  for  laser  imprint  suppression  and  adiabat 
shaping  with  the  aid  of  a  spike  pre -pulse  producing  a  shock-rarefaction  flow.  ’  ’  Exact 
and  approximate  analytical  results  of  the  linear  stability  analysis  obtained  for  the  related 
cases  of  impulsive  loading4  and  Vishniac  instability  of  a  blast  wave,35  respectively, 
indicate  feasibility  of  this  development. 
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appendix:  derivation  of  the  laplace-transormed  solutions 

When  the  shock  wave  is  driven  by  a  piston  starting  from  t  =  0 ,  we  impose  the 
boundary  condition  (1 1)  or  (12)  at  x  =  0  for  all  time.  This  boundary  condition  requires 
the  existence  of  reflected  pressure  waves  in  the  positive  x-direction,  following  the  shock 
front.  As  a  consequence,  the  Laplace  transform  for  the  shock  pressure  perturbations  has 
to  be  obtained  from  an  inhomogeneous  functional  equation.  ’  ’  Due  to  this,  the 
mathematical  structure  of  the  solution  in  a  general  case  cannot  be  presented  by  a  simple 

closed  expression  like  that  of  Eqs.  (43)  and  (68).  The  actual  formula  for  Ps  consists  of  an 
infinite  sum  of  more  or  less  complicated  functions.  Each  term  in  the  sum  appears  as  a 
consequence  of  the  multiple  reverberations  of  sonic  waves  that  take  place  between  the 
shock  surface  and  the  piston.  For  such  boundary  condition,  we  consider  the  full  problem 
of  the  “quilt”  pre-shock  density  modulation.  After  some  algebra,  using  the  method 

described  in  Ref.  16,  it  can  be  shown  that  the  pressure  function  Ps  is  the  solution  of  a 
functional  equation:16 

ps  iq)  =\  (q)+K  (q)Ps  (q + 20 ,  X  (Ai) 

where  the  Laplace-transform  variable  s  is  related  to  q  by  s  =  sinh  q  .  The  characteristic 
that  makes  the  above  equation  more  difficult  to  solve  is  the  fact  that  the  argument  of  the 
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pressure  is  “shifted”  by  the  amount  20  s  in  the  right-hand-side  of  Eq.  (Al).  This  shift  in 

the  equation  is  necessary  to  take  into  account  the  multiple  reverberations  between  the 
shock  and  the  piston,  and  the  Doppler  shift  that  occurs  after  each  reflection  at  the  shock 
surface.24  The  functions  \{q)  and  X2 (q)  for  the  double-periodic  “quilt”  pre-shock 
density  modulation,  with  rigid  piston,  free  surface  or  isolated  shock  boundary  conditions, 
are: 


\{q)  = 


a2  (q)  +  Rpa2  (q  +  20  s) 
coshg-a,(g) 


K(q)  =  R, 


al  (q  +  20  s )  +  cosh(g  +  20 s ) 
cosh^-a^g) 


where 


ai (^)  =  a10  sinh q  +  —7^ ,  a2(q)  =  -awps0  +  .  J*20  .  (A3) 

sinhg  smh  g+q0 

Here  R  is  the  reflection  coefficient  for  sonic  waves  incident  upon  the  piston  (R  =  1 , 

- 1 ,  and  0  for  the  rigid  piston,  free  surface  and  isolated  shock  boundary  conditions, 
respectively),  and  ps0  is  the  initial  pressure  perturbation  amplitude  at  the  piston  surface. 


Its  value  for  the  rigid  piston  boundary  condition  is  obtained  from  the  perturbed  Rankine- 
Hugoniot  equations  evaluated  at  t  =  0+  : 


PsO  =£a- 


m]  +  1 


(A4) 


for  the  free  surface  boundary  condition,  ps0  =  0 ,  as  follows  from  (12),  and  for  the 


isolated  shock  boundary  condition  it  is  given  by  (48). 

The  solution  to  Eq.  (Al)  is  the  sum  of  a  solution  of  the  homogeneous  equation 
plus  a  particular  solution.  As  our  gas  has  an  ideal  EOS,  it  is  not  difficult  to  see  that  the 
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only  acceptable  solution  to  the  homogeneous  equation  is  the  trivial  zero  function.  4  An 
exact  particular  solution  of  (Al)  is  easily  obtained  by  iterations:16 

ps  (q)  =  \  (q)  +  Y,\(q+ 2  «es  >17  K  (q + 2mQs )  •  (A5) 

n= 1  m= 0 

With  another  EOS,  different  from  an  ideal  gas,  the  homogeneous  equation  could 
contribute  with  a  non-zero  solution  and  therefore,  another  modes  could  be  excited  at  the 
shock  surface.24  For  the  case  of  isolated  shock  boundary  condition,  Rp  =  0  ,  and 

therefore  all  the  terms  with  shifted  arguments  q  +  20 s  describing  reverberations  in  (A2), 
(A3)  vanish,  so  that  k2  (q)  =  0 .  Hence  in  the  right-hand  side  of  (A5)  only  the  first  term 
remains,  and  this  solution  is  reduced  to  (68),  or,  in  a  particular  case  of  the  “stripes”  pre¬ 
shock  density  modulation,  to  (43).  For  other  boundary  conditions,  the  inversion  of  (A5) 
needed  to  get  the  complete  temporal  evolution  requires  the  evaluation  of  a  sufficiently 
large  number  of  terms,  similarly  to  the  case  of  the  classical  RM  instability  studied  in  Ref. 
16,  hence  the  temporal  solution  could  be  rather  complicated.  Usually,  few  terms  are 
sufficient.  We  can  calculate  the  complete  temporal  evolution  starting  from  t  =  0+  by  two 
different  procedures:  using  the  original  functional  equation  (Al)  or  its  exact  solution 

provided  by  Eq.  (A5).  If  we  expand  Ps{q )  into  a  power  series  in  cxp(-cy)  and  substitute 
it  into  Eq.  (Al),  as  suggested  for  the  first  time  by  Zaidel,"  we  could  equate  powers  of 
exp(-q)  and  get  the  complete  temporal  solution  in  the  whole  fluid  as  a  series  of  Bessel 
functions  of  integer  order,  because  the  inverse  Laplace  transform  of  the  nlh  integer  power 
of  cxp(-cy)  is  nJn(rs)/rs ,  where  rs  is  given  by  (50).  An  alternative  way  consists  of 
truncating  Eq.  (A5)  at  some  order  and  constructing  an  integral  representation  of  the 
pressure  perturbations  at  the  shock  front,  as  was  done  for  the  isolated  shock  boundary 
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condition,  see  Eqs.  (49),  (55)  and  (58).  It  must  be  taken  into  account  that  the  order  of 
truncation  will  depend  on  the  value  of  0? ,  which  is  a  monotonically  decreasing  function 

of  Mj .  As  the  shock  strength  increases,  more  terms  should  be  retained,  at  least  to 
accurately  describe  the  sonic  reverberations  and  the  vorticity/entropy  generation  behind 
the  shock  front,  when  the  shock  is  not  very  far  from  the  surface  x  =  0  .  To  do  this,  we 
must  construct  an  auxiliary  function  /  (z)  [as  in  Eqs.  (49),  (55)  and  (58)],  that  would  be 

essentially  the  imaginary  part  of  Ps  ( iz )  evaluated  at  both  sides  of  the  imaginary  axis 
along  the  segment  [-/',/] ,  with  the  multi-valued  character  of  expressions  like 

yls2  + 1  properly  taken  into  account.  It  is  quite  clear  that  more  complicated  boundary 
conditions  result  in  more  complicated  mathematical  structure  of  the  function  fp(z )  .  In 

the  case  of  an  isolated  shock,  we  have  obtained  a  very  simple  closed  form  expression  for 
this  auxiliary  kernel.  This  may  not  be  so  for  the  rigid  piston  boundary,  as  the  algebraic 
complexity  of  f  (z)  will  depend  on  the  level  of  truncation  of  the  infinite  series  given  by 

Eq.  (A5).  Nevertheless,  despite  of  these  algebraic  complications,  the  exact  solution  given 
by  Eq.  (A5)  shows  its  power  and  beauty  in  providing  us  with  important  information  on 
the  behavior  at  large  times. 

For  the  long- wavelength  regime,  <  1 ,  it  is  easy  to  see  that  we  obtain  the  same 

final  asymptotic  behavior  as  in  the  isolated  shock  case,  given  already  in  Eqs.  (32)  and 
(61).  The  only  substantial  difference  between  the  isolated  shock  and  the  rigid  piston  case 
lies  in  the  amplitudes  of  the  temporal  decaying  terms.  But  once  this  transient  behavior 
becomes  negligible,  the  final  oscillatory  response  is  the  same  for  both  cases, 
independently  of  the  boundary  conditions  far  downstream  of  the  shock  surface.  In  this 
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regime,  the  only  term  that  contributes  with  an  oscillatory  asymptotic  is  the  first  term  that 
appears  in  \(q) : 


a 


20 


1 


sinlrg+^o  coshg-a^)' 


(A6) 


with  two  poles  at  5  =  ±i£0.  This  expression  reproduces  the  exact  solution  to  the  isolated 
shock  problem,  Eq.  (68),  which  explains  why  the  asymptotic  expressions  coincide  in  the 
long-wavelength  regime.  The  other  terms  that  appear  in  Eq.  (A5)  are  due  to  the  multiple 
reverberations  that  occur  between  the  rigid  wall  and  the  shock  surface  and  their 
contribution  for  large  times  is  an  oscillation  with  an  amplitude  that  decays  as  A3 12 .  These 
additional  terms  are  actually  negligible  for  very  weak  shocks  because 
0v  =  -(l/2)ln(Mj  —  l)— >  00  as  Mx  approaches  unity,  making  those  terms  less  important 
than  the  one  shown  in  Eq.  (A6). 

However,  within  the  remaining  terms  that  represent  the  multiple  reverberations  of 
the  sonic  waves  in  Eq.  (69),  we  recognize  contributions  of  the  form: 


a 


20 


1 


sinh2 ( q  +  20  s)  + 1, l  cosh(g  +  20s) -ax(q  +  20 s) 


(A7) 


Higher-order  terms  like  one  shown  in  Eq.  (A7)  (involving  functions  with  arguments  of 
the  shifted  form  q  +  4Qs ,  q  +  60s ,  etc.)  might  contribute  pure  imaginary  poles.  These 

poles  are  imaginary  roots  of  the  equation  sinh2  (<7  +  2//0J  +(^ 2  =  0.  with  n  >  1  any 
positive  integer;  they  add  “extra”  non-decaying  oscillations  at  normalized  frequencies 
lower  than  C,0.  Indeed,  we  find  two  such  terms,  with  a  shift  equal  to  20s  for  n  =  1 .  Those 
tenns  contribute  pure  imaginary  poles  in  the  short-wavelength  regime,  when  not  only 
C,0  >  1 ,  but  also  a  stricter  condition  (77)  is  satisfied.  The  number  N  of  these  poles  is  given 
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by  Eq.  (80).  Non-decaying  oscillations  corresponding  to  the  shock  front  perturbation 
amplitudes  from  the  poles  are  given  by  the  n  -  1,2,  ...,7V  terms  of  the  sums  (81),  (82). 
The  value  of  n  for  each  of  these  terms  corresponds  to  the  number  of  reflections  from  the 
piston. 


For  the  free  surface  boundary  condition,  the  piston  surface  x  =  0  after  t  =  0+ 
becomes  a  vortex  sheet,  which  exhibits  a  RM-type  instability. 1  It  is  not  difficult  to  arrive 
to  a  formula  for  the  asymptotic  growth  rate  at  the  rippled  surface/interface,  which  is  valid 
both  for  the  isolated  shock  case  and  for  the  free  surface  case  in  the  long-wavelength 
regime.  The  calculations  are  similar  to  those  outlined  in  detail  in  Ref.  16.  We  obtain  the 
following  expression  for  the  Laplace  transform  of  the  ripple  x-velocity  at  x  =  0  : 


vxM)  = 


l  . 

- K 

sinhg 


-o+ 


[cosh(^r  +0V)  +  (X,  {q  +0S )]ps  (q  +0S )- 


a 


20 


sinh2(g+0j  +  ^( 


(A8) 


where  the  velocity  perturbation  at  t  =  0+  is  defined  below  Eq.  (88)  .  The  Laplace 
transform  (A8)  shows  a  pole  at  5  =  0,  giving  rise  to  a  constant  asymptotic  rate  of  linear 
growth. 1 6  It  also  can  have  imaginary  poles  at  L,  l ,  L,  3  . . . ,  describing  sonic  oscillations, 

when  the  conditions  for  those  oscillations  to  exist  are  satisfied:  o0  >  (2 n  - 1)0  s ,  n  —  1,2, 
. . .,  cf.  (78).  Their  number  is  given  by  (89),  and  their  nonnalized  frequencies  -  by  (90). 
The  asymptotic  time  dependence  corresponding  to  (A8)  is  given  by  (91). 
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Figure  1.  Density  maps  showing  examples  of  the  pre-shock  density  modulation: 
(a)  periodic  in  the  lateral  direction  (“stripes”);  (b)  periodic  in  both  lateral  and  longitudinal 
directions  (“quilt”);  (c)  periodic  in  the  lateral  direction  and  concentrated  near  the  surface; 
(d)  non-periodic,  Gaussian  density  profde  in  the  lateral  direction. 
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Figure  2.  Density  perturbation  map  illustrating  formulation  of  the  problem  for  the 
“stripes”  case  of  Fig.  1(a).  The  shock  wave  is  driven  by  a  rigid  piston  in  the  positive  x 
direction.  Density  modulation  of  the  pre-shock  fluid  1  is  periodic  in  y,  independent  of 
time  and  the  longitudinal  coordinate  x.  The  flow  of  the  shocked  fluid  2  is  non-steady  and 
x-dependent  but  maintains  the  original  periodicity  in  the  lateral  direction  y.  The  shape  of 
the  perturbed  shock  front  is  a  single-mode  sine  wave  whose  amplitude  is  denoted  by  8xs . 
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Figure  3.  Normalized  shock  ripple  amplitude  k8xs  vs.  normalized  time  kDt  for 
the  “stripes”  case  of  the  pre-shock  density  modulation,  y  =  5/3 ,  M{  =10 .  Dotted  line 
displays  the  asymptotic  value  (34),  thick,  thin  and  gray  solid  lines  correspond  to  the  rigid 
piston,  free  surface  and  isolated  shock  boundary  conditions,  respectively. 
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Figure  4.  Normalized  shock  displacement  amplitude  qbxs  (a)  and  normalized 
curvature  of  the  shock  front  VqR  at  r  =  0  (b)  vs.  normalized  time  qDt  for  the 
cylindrical  Gaussian  “tube”  localized  pre-shock  density  perturbation,  y=5/3,M,=10. 
Thick  and  thin  solid  lines  correspond  to  the  rigid  piston  and  free  surface  boundary 
conditions  at  the  piston  surface,  respectively.  Horizontal  dotted  lines  show  the  asymptotic 
values  (60)  and  (61). 
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Figure  5.  Normalized  shock  displacement  amplitude  qhxs  (a)  and  normalized 
curvature  of  the  shock  front  1  /  qR  at  y  =  0  (b)  vs.  normalized  time  qDt  for  the  planar 
Gaussian  “joint”  localized  pre-shock  density  perturbation,  y=5/3,M,=10.  Thick  and 
thin  solid  lines  correspond  to  the  rigid  piston  and  free  surface  boundary  conditions  at  the 
piston  surface,  respectively.  Dotted  lines  in  (a)  show  the  asymptotic  formulae  (66)  for 
each  boundary  condition;  horizontal  dotted  line  in  (b)  shows  the  asymptotic  curvature 
(67). 
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Figure  6.  Normalized  shock  ripple  amplitude  kbxs  vs.  normalized  time  kDt  for 
the  long-wavelength  “quilt”  case  of  the  pre-shock  density  modulation,  y  =5/3,  M{  =10 , 
kx  Ik  =  1/5  .  Dotted  line  displays  the  asymptotic  formula  (35),  thick,  thin  and  gray  solid 

lines  correspond  to  the  rigid  piston,  free  surface  and  isolated  shock  boundary  conditions, 
respectively. 
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Figure  7.  Normalized  shock  ripple  amplitude  kbxs\s.  normalized  time  kDt  for 
the  short-wavelength  “quilt”  case  of  the  pre-shock  density  modulation:  y  =  5  /  3 , 

My  =  10  ,  kx  /  ky  =  1 :  (a)  Dotted  line  displays  the  asymptotic  fonnula  (35),  gray  line  -  the 

exact  solution  for  isolated  shock,  (b)  Thick  and  thin  solid  lines  correspond  to  the  exact 
solution  for  rigid  piston  and  the  free  surface  boundary  conditions,  respectively,  circles 
and  diamonds  -  to  Eq.  (82)  with  R  =  1  and  - 1 ,  respectively. 
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Figure  8.  Normalized  shock  ripple  amplitude  kbxs  vs.  normalized  time  kDt  for 
the  pre-shock  standing  sonic  wave,  y  =  5/3 ,  Mx  =  10  (a)  and  Mx  =  25  (b).  Dotted  line 
displays  the  asymptotic  shape  (32),  thick  and  thin  solid  lines  correspond  to  the  rigid 
piston  and  free  surface  boundary  conditions  at  the  piston  surface,  respectively. 
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Figure  9.  Longitudinal  velocity  perturbation  amplitude  at  x  =  0  vs.  normalized 
time  kDt  for  an  isolated  shock  wave  interacting  with  a  “quilt”  pre-shock  density 
perturbation  shown  as  gray  line,  (a)  Long-wavelength  regime,  kx  / ky  =  1/5  ;  horizontal 

dashed  line  -  late-time  asymptotic  value  given  by  (88).  (b)  Short-wavelength  regime, 
kx  / k  =  1 ;  dotted  line  -  late-time  asymptotic  oscillations  given  by  (91),  horizontal 

dashed  line  -  average  late-time  value  given  by  (88).  For  both  cases,  y=5/3,M1=10. 
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Figure  10.  Normalized  displacement  amplitude  of  the  piston  surface  k8xp  (thick), 
its  time  derivative  Vp  (thin),  and  the  constant  asymptotic  value  (88)  (horizontal  dashed 

line)  (a);  time  derivative  of  areal  mass  modulation  amplitude  8m  (thick)  and  contribution 
to  it  due  to  the  surface  ripple  (thin  line)  (b)  vs.  nonnalized  time  kDt  for  the  “stripes” 
pre-shock  density  perturbation,  y=5/3,M1=10,  free  surface  boundary  conditions  at 
the  piston. 
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Figure  11.  (a)  Normalized  areal  mass  modulation  amplitude  vs.  nonnalized  time 
kDt  for  the  “stripes”  pre-shock  density  perturbation, y  =  5/3 ,  M ,  =  10  ,  rigid  piston,  (b) 
Relative  post-shock  density  modulation  amplitude  profdes  p  =  Sp  /  p2  shown  at  the 
times  when  the  shock  wave  traveled  1,  2,  3,  4  and  5  lateral  perturbation  wavelengths 
(solid  lines).  Dotted  horizontal  line  indicates  the  pre-shock  relative  density  perturbation 
amplitude. 
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Figure  12.  Normalized  displacement  of  the  ablation  front  kbxa  (a)  and  areal  mass 
modulation  amplitude  8 m  (b)  vs.  normalized  time  kDt  for  a  “stripes”  density 
perturbation  y=5/3,M,=10,  and  an  ablation  front  with  rD  =  0.0625 ,  Ma  =  0.05 
acting  as  a  piston. 
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Figure  13.  Normalized  time  derivatives  of  the  areal  mass  modulation  amplitude 
8m  for  surface  ripples  (thick  lines),  exponential  (a)  and  Bessel-function  (b)  localized 
density  perturbations  in  the  target  (thin  lines,  labeled  by  the  value  of  the  parameter  q/k) 
vs.  nonnalized  time  kDt  for  y  =  5/3 ,  Mx  =  10  . 
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Figure  14.  Theoretical  and  simulated  time  histories  of  normalized  areal  mass 
modulation  amplitude  for  the  “stripes”  case:  y  =  5/3 ,  M ,  =  10  ,  lateral  perturbation 

wavelength  X  =  30  pm,  rigid  piston  velocity  V  =  4xl06  cm/s.  Lines  1  through  4 

correspond  to  the  values  of  the  small-amplitude  parameter  £  =  5xl0^4 ,  2.5x10  5 , 

1 0  2  and  5xl()  2 ,  respectively,  (a)  Simulation  done  with  a  constant  time  step;  (b)  time 
step  initially  reduced. 
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Figure  15.  Normalized  areal  mass  modulation  amplitude  8m  vs.  time  for  density 
modulation  profiles  (97)  localized  near  the  surface:  y  =  5/3 ,  Mx  =  10 ,  lateral 

perturbation  wavelength  X  =  2n  / q  =  30  pm,  rigid  piston  velocity  V  =  4x  106  cm/s;  in 

the  simulations,  effective  ripple  amplitude  8x0  =  0.1  pm:  (a)  theory,  (b)  simulation  with 

constant  time  step,  (c)  simulations  with  initially  reduced  time  step.  Lines  0  through  4  are 
indexed  by  the  value  of  the  parameter  n  in  (97). 
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